Decomposing Treatment Effect Variation
Peng Ding
UC Berkeley

Avi Feller
UC Berkeley

∗

Luke Miratrix
Harvard GSE

arXiv:1605.06566v2 [math.ST] 28 Jul 2017

July 31, 2017

Abstract
Understanding and characterizing treatment effect variation in randomized experiments has
become essential for going beyond the “black box” of the average treatment effect. Nonetheless,
traditional statistical approaches often ignore or assume away such variation. In the context of
randomized experiments, this paper proposes a framework for decomposing overall treatment
effect variation into a systematic component explained by observed covariates and a remaining
idiosyncratic component. Our framework is fully randomization-based, with estimates of treatment effect variation that are entirely justified by the randomization itself. Our framework can
also account for noncompliance, which is an important practical complication. We make several contributions. First, we show that randomization-based estimates of systematic variation
are very similar in form to estimates from fully-interacted linear regression and two stage least
squares. Second, we use these estimators to develop an omnibus test for systematic treatment
effect variation, both with and without noncompliance. Third, we propose an R2 -like measure
of treatment effect variation explained by covariates and, when applicable, noncompliance. Finally, we assess these methods via simulation studies and apply them to the Head Start Impact
Study, a large-scale randomized experiment.
Key Words: Noncompliance; Heterogeneous treatment effect; Idiosyncratic treatment effect
variation; Randomization inference; Systematic treatment effect variation.

∗

Peng Ding (Email: pengdingpku@berkeley.edu) is Assistant Professor, Department of Statistics, University of
California, Berkeley. Avi Feller (Email: afeller@berkeley.edu) is Assistant Professor, Goldman School of Public
Policy, University of California, Berkeley. Luke Miratrix (Email: lmiratrix@g.harvard.edu) is Assistant Professor,
Harvard Graduate School of Education. We thank Alberto Abadie, Donald Rubin, participants at the Applied
Statistics Seminar at the Harvard Institute of Quantitative Social Science, and colleagues at University of California,
Berkeley and Harvard University for helpful comments. We gratefully acknowledge financial support from the Spencer
Foundation through a grant entitled “Using Emerging Methods with Existing Data from Multi-site Trials to Learn
About and From Variation in Educational Program Effects,” and from the Institute for Education Science (IES Grant
#R305D150040).

1

1

Introduction

The analysis of randomized experiments has traditionally focused on the average treatment effect,
often ignoring or assuming away treatment effect variation (e.g., Neyman, 1923; Fisher, 1935;
Kempthorne, 1952; Rosenbaum, 2002). Today, understanding and characterizing treatment effect
variation in randomized experiments has become essential for going beyond the “black box” of
the average treatment effect. This is clear from the increasing number of papers on the topic in
statistics and machine learning (Hill, 2011; Athey and Imbens, 2016; Wager and Athey, 2017),
biostatistics (Huang et al., 2012; Matsouaka et al., 2014), education (Raudenbush and Bloom,
2015), economics (Heckman et al., 1997; Crump et al., 2008; Djebbari and Smith, 2008), political
science (Green and Kern, 2012; Imai and Ratkovic, 2013), and other areas.
This paper proposes a framework for decomposing overall treatment effect variation in a randomized experiment into a systematic component that is explained by observed covariates, and an
idiosyncratic component that is not explained (Heckman et al., 1997; Djebbari and Smith, 2008). In
doing so, we make several key contributions. First, we take a fully randomization-based perspective
(see Rosenbaum, 2002; Imbens and Rubin, 2015), and propose estimators that are entirely justified
by the randomization itself. This is in contrast to much of the literature on randomization-based
methods, where treatment effect variation is typically a nuisance (e.g. Rosenbaum, 1999, 2007).
Similar to Lin (2013), we show that the resulting estimator is very similar in form to linear regression with interactions between the treatment indicator and covariates. Unlike with linear regression,
however, the proposed estimator does not require any modeling assumptions on the marginal outcomes.
Second, we extend these methods from intention-to-treat (ITT) analysis to allow for noncompliance, proposing a randomized-based estimator for systematic treatment effect variation for the
Local Average Treatment Effect (LATE) in the case of noncompliance (Angrist et al., 1996). We
show that this estimator is nearly identical to the two-stage least squares estimator with interactions between the treatment and covariates. We believe that this is a particularly novel contribution
to the recent literature seeking to reconcile the randomization-based tradition in statistics and the
linear model-based perspective more common in econometrics (Abadie, 2003; Imbens, 2014; Imbens
and Rubin, 2015).
Armed with these estimators, we turn to two practical tools for decomposing treatment effect
2

variation. The first is an omnibus test for the presence of systematic treatment effect variation.
While versions of this test have been proposed previously, largely in the context of linear models (Cox, 1984; Crump et al., 2008), our proposed test is fully randomization-based and can also
account for noncompliance. The second is to develop and bound an R2 -like measure of the fraction
of treatment effect variation explained by covariates. This builds on previous versions proposed
in the econometrics literature (Heckman et al., 1997; Djebbari and Smith, 2008), again extending
results to account for noncompliance. This approach is also closely related to the Oaxaca–Blinder
decomposition in economics (Oaxaca, 1973; Blinder, 1973). See Angrist et al. (2013) for a recent
application that also addresses compliance. Finally, we apply these methods to the Head Start
Impact Study, a large-scale randomized trial of Head Start, a federally funded preschool program
(Puma et al., 2010). We relegate the technical details and some further extensions to the online
Supplementary Material.

2

Framework for Treatment Effect Variation

2.1

Setup and notation

Assume that we have n units in an experiment. For unit i, let Xi = (X1i , . . . , XKi )T ∈ RK denote
the vector of pretreatment covariates, with the constant 1 as its first component. Let Ti denote
the treatment indicator with 1 for treatment and 0 for control. We use the potential outcomes
framework (Neyman, 1923; Rubin, 1974) to define causal effects. Under the Stable Unit Treatment
Value Assumption (Rubin, 1980) that there is only one version of the treatment and no interference
among units, we define Yi (1) and Yi (0) as the potential outcomes of unit i under treatment and
control, respectively. The observed outcome, Yiobs = Ti Yi (1) + (1 − Ti )Yi (0), is quite general
and includes continuous, binary, and zero-inflated cases. On the difference scale, the individual
treatment effect is τi = Yi (1) − Yi (0).
Importantly, this is finite population inference in that we condition on the n units at hand—the
potential outcomes are fixed and pre-treatment. This differs from super population inference in
which some variables or residuals are assumed to be independent and identically distributed (iid)
draws from some distribution. See, for example, Rosenbaum (2002), Imbens and Rubin (2015)
and Li and Ding (2017). Under the potential outcomes framework, {Yi (1), Yi (0)}ni=1 are all fixed
numbers; the randomness of any estimator comes from the assignment mechanism, which is the

3

distribution of possible treatment assignments T = (T1 , . . . , Tn )T . Note that pr{(T1 , . . . , Tn ) =
−1 Pn
(t1 , . . . , tn )} = nn1
if i=1 ti = n1 .

2.2

Randomization inference for vector outcomes

To set up our overall framework, we first generalize Neyman (1923)’s classic results to vector
outcomes. We consider a completely randomized experiment, with n1 units assigned to treatment

and n0 units assigned to control; in total we have nn1 possible randomizations. We are interested
in estimating the finite population average treatment effect on a vector outcome V ∈ RK :
n

τV =

1X
{Vi (1) − Vi (0)} ,
n
i=1

where Vi (1) and Vi (0) are the potential outcomes of V for unit i. The Neyman-type unbiased
estimator for τV is the difference between the sample mean vectors of the observed outcomes under
treatment and control:
τbV =

V̄1obs

−

V̄0obs

n
n
n
n
1 X
1 X
1 X
1 X
obs
obs
=
Ti Vi −
(1 − Ti )Vi =
Ti Vi (1) −
(1 − Ti )Vi (0).
n1
n0
n1
n0
i=1

i=1

i=1

i=1

The behavior of our estimator, and of our estimators for heterogeneity discussed later, revolve
around covariances of vector outcomes. For notation, let A = {A1 , . . . , An } be a collection of n
P
vectors, with Ā = n−1 ni=1 Ai the vector mean, and define the covariance operator on A as
n

S(A) =

1 X
(Ai − Ā)(Ai − Ā)T ,
n−1
i=1

which gives the covariance matrix of the n vectors in A. For example, Ai can be Vi (1), Vi (0) or
Vi (1) − Vi (0).
The following theorem, generalizing the results for scalar outcomes from Neyman (1923), demonstrates that τbV is unbiased and gives its covariance matrix.
Theorem 1. Over all possible randomizations of a completely randomized experiment, τbV is
unbiased for τV , with K × K covariance matrix:
cov(b
τV ) =

S{V (1)} S{V (0)} S{V (1) − V (0)}
+
−
.
n1
n0
n

(1)

The diagonal elements of this matrix are the variances of the estimators of each component of
τV . The covariance matrix of τbV depends on the various covariances of the potential outcomes
4

under treatment and control. In particular, the last term depends on the correlation between the
potential outcomes V (1) and V (0), and therefore cannot be identified from the observed data.
When the individual treatment effects are constant for all components of V , the last term in the
above covariance matrix vanishes, because then S{V (1) − V (0)} = 0K×K . Under this assumption,
we can unbiasedly estimate the sampling covariance matrix cov(b
τV ) by replacing the covariances
of the potential outcomes by the sample analogues:
cov(b
c τV ) =

Sb1 (V obs ) Sb0 (V obs )
+
,
n1
n0

where
n

Sbt (V obs ) =

1 X
I(Ti =t) (Vi − V̄tobs )(Vi − V̄tobs )T
nt − 1

(t = 0, 1)

(2)

i=1

are the sample covariance matrices of V obs in the treatment and control groups. Without the
constant treatment effect assumption, the covariance estimator cov(b
c τV ) is conservative in the sense
that the difference between the expectation of the variance estimator and the true variance is a
non-negative definite matrix. In particular, the diagonal terms of the expected estimator will all
be larger than the truth. Letting K = 1, the covariance matrices become simple variances, which
recovers Neyman’s original result.
Using the mathematical framework introduced in the Appendix and in Li and Ding (2017), we
can easily generalize Theorem 1 to more complicated experimental designs, e.g., cluster-randomized
trials (Middleton and Aronow, 2015) and unbalanced 22 split-plot designs (Zhao et al., 2017).

2.3

Decomposing Treatment Effect Variation

We now apply this general framework to treatment effect variation. We decompose the individual
treatment effect, τi , via
τi = Yi (1) − Yi (0) = XiT β + εi ,

(i = 1, . . . , n)

(3)

with β being the finite population linear regression coefficient of τi on Xi , defined by
β = arg min

b∈RK

n
X

(τi − XiT b)2 .

(4)

i=1

Following Heckman et al. (1997) and Djebbari and Smith (2008), we call δi = XiT β the systematic
treatment effect variation explained by the observed covariates, Xi , and call εi the idiosyncratic
treatment effect variation not explained by Xi .
5

More generally, we can view this decomposition in a regression-style framework. Define
n

Sxx =

1X
Xi XiT ∈ RK×K ,
n

n

Sxε =

i=1

1X
εi Xi ∈ RK ,
n

n

Sxτ =

i=1

1X
τi Xi ∈ RK ,
n
i=1

where Sxx is non-degenerate, analogous to the usual full rank assumption in linear models. Also
define

n

Sxt

1X
=
Xi Yi (t) ∈ RK ,
n

(t = 0, 1).

i=1

These are all finite population quantities, as in they are fixed pre-randomization values. The
definition of β gives Sxε = 0, i.e., εi and Xi have finite population covariance zero. Therefore, in the
spirit of the agnostic regression framework (e.g., Lin, 2013), the systematic component, δi = XiT β,
is a projection of τi onto the linear space spanned by Xi , and the idiosyncratic treatment effect,
εi , is the corresponding residual. The linear projection applies to general outcomes, including the
binary case.
Because of our finite population focus, if we observed all the potential outcomes we could
immediately calculate all individual treatment effects and apply standard linear regression theory
to (3) and obtain β. In particular, the solution of (4), i.e. the ordinary least squares (OLS) solution
from regressing τ on X, is
−1
−1
−1
β = Sxx
Sxτ = Sxx
Sx1 − Sxx
Sx0 ≡ γ1 − γ0 ,

(5)

−1 S
−1
where γ1 = Sxx
x1 and γ0 = Sxx Sx0 are the corresponding finite population regression coefficients

of the potential outcomes on the covariates. Let ei (1) = Yi (1) − XiT γ1 and ei (0) = Yi (0) − XiT γ0 be
the residual potential outcomes from the regression of Yi (t) onto X. Our idiosyncratic treatment
variation is then the difference of residuals: εi = ei (1) − ei (0). In practice, we do not fully observe
these components, but we can obtain unbiased or consistent estimates for them as we discuss below.

3
3.1

Systematic treatment effect variation for the ITT
Randomization-based estimator

We now turn to estimating β. As shown in (5), β has three components. The first term, Sxx ,
is fully observed as all the covariates are observed. Our estimation then depends on the sample
analogues of Sx1 and Sx0 :
n
1 X
b
Ti Yiobs Xi ∈ RK ,
Sx1 =
n1
i=1

n
1 X
b
Sx0 =
(1 − Ti )Yiobs Xi ∈ RK .
n0
i=1

6

bxt ’s capture how the observed potential outcomes correlate with the covariates. Plug these
The S
into (5) to obtain an overall estimate of β. The randomization of T then justifies the following
theorem.
−1 S
bx1 and S −1 S
b
Theorem 2. Under decomposition (3), Sxx
xx x0 are unbiased estimates of γ1 and γ0 ,

respectively. Therefore
−1 b
−1 b
βbRI = Sxx
Sx1 − Sxx
Sx0 ,

is an unbiased estimator for β with covariance matrix


S{Y (0)X} S(τ X)
−1
−1 S{Y (1)X}
b
+
−
Sxx
.
cov(βRI ) = Sxx
n1
n0
n

(6)

Here, for example, S{Y (0)X} denotes the covariance operator on new unit-level variables
Yi (0)Xi ∈ RK , made by scaling the Xi vector of each unit by Yi (0), similarly for S{Y (1)X}
and S(τ X). This slight abuse of notation gives formulae less cluttered by subscripts and excessive
annotation. As with the vector version of Neyman’s formula, the square root of the diagonal of
cov(βbRI ) gives the standard errors of β.
The covariance formula (6) generalizes the result of Neyman (1923) for the average treatment
effect, reducing to Neyman’s formula if Xi = 1 for all units. We can obtain a “conservative”
estimate of cov(βbRI ) by
"
cov(
c βbRI ) =

−1
Sxx

#
Sb1 (Y obs X) Sb0 (Y obs X)
−1
+
Sxx
,
n1
n0

recalling the definitions of the sample covariance operators Sb1 and Sb0 introduced in (2). Similar
to Neyman (1923), this implicitly assumes S(τ X) = 0. Under the assumption that εi = 0 for all
units (i.e., no idiosyncratic variation whatsoever), we can instead use S(b
τ X) with τb = XiT βbRI as a
plug-in estimate for S(τ X). This yields tighter standard errors based on the diagonal elements of
the covariance matrix.
Finite population asymptotic analysis

Theorem 2 holds for any finite sample. To obtain

confidence intervals and to conduct hypothesis testing as we describe below, we need to prove further that βbRI is asymptotically normal with mean β and covariance cov(βbRI ). Finite population
asymptotic analysis, however, has a slightly different flavor from the usual super population approach. Formally, the finite asymptotic scheme embeds the finite population {(Xi , Yi (1), Yi (0))}ni=1
7

with size n into a hypothetical sequence of finite populations with sizes approaching infinity. This
effectively assumes that all the finite population quantities, for example, Sxx and β, depend on
n, although they are fixed numbers for a given finite population. Moreover, the sample quantities
bx1 and βbRI depend on n as well, and are random quantities due to the randomization of
such as S
T . For notational simplicity, we drop the index n for all these quantities. Importantly, we must
impose some regularity conditions on the hypothetical sequence of finite populations. Throughout
the paper, we invoke the following conditions for asymptotic analysis, which are required for a form
of the finite population central limit theorem discussed in Li and Ding (2017, Theorem 5).
Condition 1. (i) Stable treatment proportions: p1 = n1 /n and p0 = n0 /n have positive limiting
values; (ii) Stable means, variances and covariances: the finite population means, variances and
covariances of the covariates and potential outcomes have finite limiting values; (iii) max1≤i≤n ||Vi −
V̄ ||22 /n → 0, where Vi can be the covariate vector, the outcome, and the products of them.
Under these conditions, we can extend Theorem 2 to a sequence of finite populations:



√ 
d
n βbRI − β → N 0, lim [p1 S{Y (1)X} + p0 S{Y (0)X} − S(τ X)] .
n→∞

(7)

As a result, we can state that βbRI is approximately normal with mean β and covariance matrix
(6), which allows us to construct confidence intervals and hypothesis tests. In our theory below,
we use this informal statement instead of (7) to avoid notational complexity.
Conditions (i) and (ii) are natural. Condition (iii) holds if V has more than two moments (Li
and Ding, 2017). For bounded covariates and outcomes, (iii) is satisfied automatically. For more
technical discussion of finite population causal inference, see Ding (2014), Aronow et al. (2014),
and Middleton and Aronow (2015); for regularity conditions of the finite population central limit
theorems, see Hájek (1960) and Lehmann (1998). A recent review is Li and Ding (2017).

3.2

Regression with treatment-covariate interactions

The results from randomization inference can shed light on the familiar case of linear regression
with treatment-covariate interactions. This classical approach assumes the model
Yiobs = XiT γ + Ti XiT β + ui ,

(i = 1, . . . , n)

(8)

where {ui }ni=1 are errors implicitly assumed to induce the randomness, and where β models systematic treatment effect variation, as in (3). Departing from much of the previous literature (e.g.,
8

Cox, 1984; Berrington de González and Cox, 2007; Crump et al., 2008), we study the properties
of the least squares estimator under complete randomization, without assuming that model (8) is
correctly specified. In particular, we do not assume any i.i.d. sampling; the assignment mechanism
drives the distribution of the OLS estimator.
Theorem 3. The OLS estimator for β from fitting model (8) can be rewritten as
b−1 S
b
b−1 b
βbOLS = S
xx,1 x1 − Sxx,0 Sx0 ,
where

n
X
bxx,t = 1
S
I(Ti =t) Xi XiT ,
nt

(t = 0, 1).

i=1

b−1 S
b
b−1 b
Over all possible randomizations of T , S
xx,1 x1 and Sxx,0 Sx0 are consistent estimates of γ1 and γ0
respectively; βbOLS therefore follows an asymptotic normal distribution with mean β and covariance
matrix
−1
cov(βbOLS ) = Sxx




S{e(1)X} S{e(0)X} S(εX)
−1
+
−
Sxx
.
n1
n0
n

(9)

with ei (1), ei (0), and εi as defined after (5).
b−1 S
b
b−1 S
b
b1,OLS = S
b0,OLS = S
This estimate is simply the difference between γ
xx,1 x1 and γ
xx,0 x0 ,
two OLS regressions run separately on each treatment arm. For treated units, define residual
b1,OLS , and for control units, define residual ebi = Yiobs − XiT γ
b0,OLS . We can drop
ebi = Yiobs − XiT γ
the unidentifiable term S(εX), estimate S{e(1)X} and S{e(0)X} by their sample analogues, and
conservatively estimate the asymptotic covariance matrix (9) by
"
"
#
#
b1 (b
b0 (b
S
e
X)
S
e
X)
−1
−1
−1
b
b
b
b−1 .
cov(
c βbOLS ) = S
S
S
xx,1
xx,1 + Sxx,0
xx,0
n1
n0
This form of the sandwich variance estimator has the same probability limit as the HuberWhite
covariance estimator for linear model (8) (Huber, 1967; White, 1980; Lin, 2013; Angrist and Pischke,
2008).
Importantly, βbRI and βbOLS are quite similar in form. In particular, βbRI uses the true Sxx while
bxx,0 and S
bxx,1 . The latter
βbOLS separately estimates the covariance matrix for each treatment arm, S
is effectively a ratio estimator. Although this introduces some small bias (on the order of 1/n),
bxx,t rather than true Sxx can often lead to gains in precision, especially when
using the estimated S
9

covariates are strongly correlated with the potential outcomes. In particular, the OLS estimator,
by separately estimating the (known) Sxx matrix for each treatment arm, can account for random
imbalances in the covariates in both arms.
The RI estimator, by comparison, has no adjustment whatsoever, and so cannot account for such
random covariate imbalances. However, in Section 3.4 below, and in the supplementary materials
we introduce a different form of adjustment that uses covariates to make the estimates of the Sxt
more precise. Depending on the structure of covariates, this estimator could be better or worse
than OLS adjustment; we leave a thorough investigation of these trade-offs for future work.
Regardless, we again emphasize that we do not rely on classical OLS assumptions to justify
the OLS estimator here. Rather, randomization (plus some mild regularity conditions for the
finite sample asymptotics) justify our results. For related discussion, see Cochran (1977) on ratio
estimators in surveys.

3.3

Omnibus test for systematic variation

Finally, we can use these results to develop an omnibus test for the presence of any systematic
treatment effect variation. The null hypothesis of no treatment effect variation explained by the
observed covariates can be characterized by
H0 (X) : β1 = 0,
where β1 contains all the components of β except the first component corresponding to the intercept. Under H0 (X), the individual treatment effects have no linear dependance on X.
We then construct a Wald-type test for H0 (X) using an estimator βb and its covariance estimator
b it could be βbRI or βbOLS . Let βb1 and cov(
cov(
c β);
c βb1 ) denote the sub-vector of βb and sub-matrix of
b corresponding to the non-intercept coordinates of X. We reject when
cov(
c β),
βb1T cov
c −1 (βb1 )βb1 > qK−1 (1 − α),

(10)

where qK−1 (1 − α) is the 1 − α quantile of the χ2 random variable with degrees of freedom K − 1.
The test in (10) is nearly identical to the test proposed by Crump et al. (2008). They relax the
parametric assumption by taking a “sieve estimator” approach, namely by using a quadratic form of
the regression function, which allows for more flexible marginal distributions. Our approach differs
in that we avoid modeling the marginal distributions entirely. If desired, we can add polynomials
10

of X (or other basis functions) into the model for δ to allow for more flexible systematic treatment
effect variation, which could enhance power or model more complex relationships between the bX
and treatment impact.

3.4

Additional considerations

In the Supplementary Material, we describe two additional points about systematic treatment
effect variation that we briefly address here. First, as mentioned above, we can use model-assisted
estimation to improve the randomization-based estimator. In particular, improving estimation
bxt directly improves βbRI , as the S
bxt are the only random components. In particular, if we
of S
bxt , by a more efficient, model-assisted estimator, as in
replace the standard sample estimator, S
survey sampling (Cochran, 1977; Särndal et al., 2003), we can achieve meaningful precision gains
in practice. More importantly, this setup allows researchers to assess systematic variation across
one set of covariates while adjusting for another set.
Second, under the assumption of no idiosyncratic variation (i.e., εi = 0 for all i), we can obtain exact inference for β by inverting a sequence of randomization-based tests. This complements
previous work on randomization-based tests for the presence of idiosyncratic treatment effect variation (Ding et al., 2016).

4

Idiosyncratic treatment effect variation for ITT

After characterizing the systematic component of treatment effect variation, we now turn to characterizing the idiosyncratic component. Since this quantity is inherently unidentifiable, we propose
sharp bounds on this component and a framework for sensitivity analysis. We then leverage these
results to bound an R2 -like measure of the treatment effect variation explained by covariates.

4.1

Bounds

We first define the main quantities of interest:
n

Sτ τ =

1X
(τi − τ )2 ,
n
i=1

n

Sδδ =

1X
(δi − τ )2 ,
n
i=1

n

Sεε =

1X 2
εi ,
n
i=1

with δi and εi defined as in (3). Then Sτ τ = Sδδ + Sεε . We can immediately estimate Sδδ via
b n , where βb is a consistent estimator, e.g., βbRI or βbOLS .
the sample variance of {δbi = XiT β}
i=1

11

However, the idiosyncratic variance, Sεε , is inherently unidentifiable because it depends on the
joint distribution of potential outcomes.
We can, however, derive sharp bounds for Sεε . Let F1 (y) and F0 (y) be the empirical cumulative
distribution functions of {ei (1)}ni=1 and {ei (0)}ni=1 . Below we denote e(t) as a random variable taking equal probabilities on n values of {ei (t)}ni=1 . Based on the Fréchet–Hoeffding bounds (Hoeffding,
1941; Fréchet, 1951; Nelsen, 2007), we can bound Sεε as follows.
Theorem 4. Sεε has sharp bounds S εε ≤ Sεε ≤ S εε , where
Z 1
Z 1
−1
−1
2
{F1 (u) − F0 (u)} du, S εε =
{F1−1 (u) − F0−1 (1 − u)}2 du
S εε =
0

with

F −1 (u)

0

= inf{x : F (x) ≥ u} as the quantile function. The lower and upper bounds are

attainable when e(1) and e(0) have the same ranks and opposite ranks, respectively.
The lower bound of Sεε corresponds to a rank-preserving relationship between e(1) and e(0),
and the upper bound of Sεε corresponds to an anti-rank-preserving relationship between e(1) and
e(0). Equivalently, they correspond to the cases where the Spearman rank correlation coefficients
between e(1) and e(0) are +1 and −1.
In practice, we can often sharpen these bounds because we are unlikely to have negatively
associated potential outcomes after adjusting for covariates. If we assume a nonnegative correlation
between e(1) and e(0), we have the following corollary:
Corollary 1. If the correlation between e(1) and e(0) is nonnegative, then the bounds for Sεε
become S εε ≤ Sεε ≤ V1 + V0 , where Vt is the variance of e(t) for t = 0, 1.
b Fe1 (y) and
We can consistently estimate each quantity: Sδδ by the sample variance of XiT β,
Fe0 (y) by Fb1 (y) and Fb0 (y), the empirical cumulative distribution functions of the residuals ebi under
treatment and control, and V1 and V0 by the variances of eb(1) and eb(0).
Variance of the overall ITT estimator.

We can use these results to obtain sharper bounds
Pn
Pn
obs − n−1
on the variance of Neyman (1923)’s estimate of overall ITT, τb = n−1
1
0
i=1 Ti Yi
i=1 (1 −
Ti )Yiobs , extending previous work by Heckman et al. (1997) and Aronow et al. (2014). See also Fogarty (2016). Applying the results in Section 2 for scalar outcomes, we have the following variance
for the difference-in-means estimator,
S11 S00
var(b
τ) =
+
−
n1
n0
12



Sδδ
Sεε
+
n
n


,

where Sτ τ = Sδδ + Sεε . As we discuss above, Neyman (1923) proposed a lower bound for the overall
var(b
τ ) under the assumption of a constant treatment effect, Sτ τ = 0. More recently, Aronow et al.
(2014) instead proposed to bound Sτ τ via Fréchet–Hoeffding bounds. We can modestly improve
these results by applying Fréchet–Hoeffding bounds for Sεε alone rather than for Sτ τ = Sδδ + Sεε .
So long as Sδδ > 0, this yields strictly tighter bounds on var(b
τ ) than the corresponding bounds
that do not incorporate covariate information. In turn, this gives a tighter estimate of the standard
error for the same difference-in-means estimator, τb.
A variance ratio test. Finally, while the relationship between e(0) and e(1) is inherently unidentifiable, there is some information in the data about the relationship between εi , the individual-level
idiosyncratic treatment effect, and Yi (0), the control potential outcome. In particular, Raudenbush
and Bloom (2015) noted that if the variance of the treatment potential outcomes is smaller than
the variance of the control potential outcomes, then the treatment effect must be negatively associated with the control potential outcomes. In the Supplementary Material, we extend this result
to incorporate covariates and propose a formal test.

4.2

Sensitivity analysis

Going beyond worst-case bounds, we can assess the sensitivity of our estimate of Sεε to different assumptions of the dependence between potential outcomes. Using the probability integral
transformation, we represent the residual potential outcomes as
e(1) = F1−1 (U1 ),

e(0) = F0−1 (U0 ),

U1 , U0 ∼ Uniform(0, 1),

Therefore, the dependence of the potential outcomes is determined by the dependence of the uniform
random variables U1 and U0 , which are the standardized ranks of the potential outcomes. When
U1 = U0 , Sεε attains the lower bound S εε ; when U1 = 1 − U0 , Sεε attains the upper bound S εε ;
when U1 U0 , Sεε attains the improved upper bound V1 + V0 .
Rather than simply examine extreme scenarios of Sεε , we can instead represent U1 as a mixture
of U0 and another independent uniform random variable V0 :
U1 ∼ ρU0 + (1 − ρ)V0 ,

i.i.d.

U0 , V0 ∼ Uniform(0, 1),

(11)

which the sensitivity parameter ρ captures the association between U1 and U0 . An immediate
interpretation of ρ is the proportion of rank preserved units, with the other 1 − ρ as the proportion
13

of units with independent treatment and control residual outcomes. When ρ = 0, U1 U0 , and
the residual potential outcomes are independent; when ρ = 1, U1 = U0 , and the residual potential
outcomes have the same ranks. The values between (0, 1) corresponds to positive rank correlation
but not full rank preservation. Note that the representation of the joint distribution is not unique,
because we can choose any copula as a joint distribution of (U1 , U0 ) (Nelsen, 2007). We choose the
above representation and notation ρ for the following theorem.
Theorem 5. If Equation 11 holds, then ρ is Spearman’s rank correlation coefficient between e(1)
and e(0). Furthermore, Sεε is a linear function of ρ:
Sεε (ρ) = ρS εε + (1 − ρ)(V1 + V0 ).
We cannot extract any information about ρ from the data. We therefore treat ρ as a sensitivity
parameter, choose a plausible range of ρ, and obtain corresponding values for Sεε .

4.3

Fraction of treatment effect variation explained

A natural question is the relative magnitudes of Sδδ and Sεε (Djebbari and Smith, 2008). Continuing
the regression analogy, this is an R2 -like measure for the proportion of total treatment effect
variation explained by the systematic component:
Rτ2 =

Sδδ
Sδδ
=
,
Sτ τ
Sδδ + Sεε

which is the ratio between the finite population variances of δ and τ. As above, we can directly
estimate Sδδ but must bound Sεε . Applying Theorem 4, we obtain the following bounds on Rτ2 .
Corollary 2. The sharp bounds on Rτ2 are
Sδδ
Sδδ
.
≤ Rτ2 ≤
Sδδ + S εε
Sδδ + S εε
If we further assume that the correlation between e(1) and e(0) is nonnegative, the sharp bounds
on Rτ2 are
Sδδ
Sδδ
≤ Rτ2 ≤
.
Sδδ + V1 + V0
Sδδ + S εε
We estimate these bounds via plug-in estimates. Note that Djebbari and Smith (2008) explore
a similar quantity by using a permutation approach to approximate the Fréchet–Hoeffding upper
and lower bounds. Finally, we can use the sensitivity results for Sεε , with values of ρ ∈ [0, 1]:
Rτ2 (ρ) =

Sδδ
.
Sδδ + Sεε (ρ)
14

5

Noncompliance

5.1

Setup

We now extend our results to allow for noncompliance. Let T be the indicator of treatment
assigned, D be the indicator of treatment received, Y be outcome of interest, and X be pretreatment
covariates. Under the Stable Unit Treatment Value Assumption, we define Di (t) and Yi (t) as
the potential outcomes for unit i under treatment assignment t. Following Angrist et al. (1996)
and Frangakis and Rubin (2002), we can classify units into four compliance types based on the
joint values of Di (1) and Di (0):


Always Taker (a)



Never Taker (n)
Ui =

Complier (c)



Defier (d)

if
if
if
if

Di (1) = 1, Di (0) = 1,
Di (1) = 0, Di (0) = 0,
Di (1) = 1, Di (0) = 0,
Di (1) = 0, Di (0) = 1.

Denote ni and πu by the number and proportion of compliance types πu of stratum U = u for
u = a, n, c, d.
Throughout our discussion, we invoke the following assumptions which are commonly used for
analyzing randomized experiments with noncompliance.
Assumption 1. (i) Monotonicity: Di (1) ≥ Di (0); (ii) Exclusion restrictions for Always Takers and
Never Takers: Yi (1) = Yi (0) for all units with Di (1) = Di (0); (iii) Strong instrument: πc > C0 > 0,
where C0 is a positive constant independent of the sample size.
Monotonicity rules out the existence of Defiers, i.e., πd = 0. Under monotonicity, we can
estimate the proportion πu using the observed counts of units classified by T and D: let ntd = #{i :
Ti = t, Di = d}, and then π
bn = n10 /n1 , π
ba = n01 /n0 , and π
bc = n11 /n1 − n01 /n0 . The exclusion
restrictions assume that treatment assignment has no effect on the outcome for Always Takers and
Never Takers. As a result, treatment effect variation is trivially zero for Always Takers and Never
Takers. Note that this is the unit-level exclusion restriction imposed in Angrist et al. (1996). This
can be relaxed in other settings; for example, we could assume the impact of randomization for these
groups is zero on average (see Imbens and Rubin, 2015). Finally, to avoid technical complexity, we
rule out the weak instrument case (Bound et al., 1995; Staiger and Stock, 1997), i.e., πc is within
a small neighborhood of 0 with radius shrinking to 0.

15

We are interested in treatment effect variation among Compliers, which motivates the following
decomposition:

τi = Yi (1) − Yi (0) =

0,
if Ui = a or n,
T
Xi βc + εi , if Ui = c,

(12)

where βc is the regression coefficient of τi on Xi among Compliers, analogous to (3).

5.2
5.2.1

Systematic treatment effect variation among Compliers
Randomization inference

We now extend the results of Section 3 to estimate systematic treatment effect variation among
Compliers. Define
Sxx,u =

n
1 X
I(Ui =u) Xi XiT ,
nu

Sxt,u =

i=1

n
1 X
I(Ui =u) Yi (t)Xi ,
nu

(t = 0, 1; u = a, c, n).

i=1

Then, analogous to (5),
−1
−1
−1
βc = Sxx,c
(Sx1,c − Sx0,c ) = Sxx,c
Sx1,c − Sxx,c
Sx0,c ≡ γ1c − γ0c ,

(13)

where
−1
γ1c = Sxx,c
Sx1,c ,

−1
γ0c = Sxx,c
Sx0,c

are the linear regression coefficients of Y (1) and Y (0) on covariates among Compliers.
Unlike in the ITT case, we cannot estimate these quantities directly. Instead, following standard
results from noncompliance (e.g., Angrist et al., 1996; Abadie, 2003; Angrist and Pischke, 2008),
we use estimates from observed subgroups to estimate the desired quantities of interest. Define
sample moments:
n
X
bxx,td = 1
I(Ti =t) I(Di =d) Xi XiT ,
S
nt
i=1

n
X
bxt,td = 1
S
I(Ti =t) I(Di =d) Yiobs Xi
nt

(t, d = 0, 1).

(14)

i=1

The following theorem connects these quantities with the finite population quantities in (13).
Theorem 6. Over all possible randomizations of a completely randomized experiment, both
bxx (1) = S
bxx,11 − S
bxx,01 and S
bxx (0) = S
bxx,00 − S
bxx,10 are unbiased for πc Sxx,c , and
S
bx1,11 − S
bx0,01 ) = πc Sx1,c ,
E(S

bx0,00 − S
bx1,10 ) = πc Sx0,c .
E(S

16

(15)

This theorem shows that we can obtain unbiased estimates for all terms in (13). The following
corollary shows that we can then obtain consistent estimates for γ1c , γ0c , and βc , recalling that in
the asymptotic analysis, we need to embed {Xi , Yi (1), Yi (0), Di (1), Di (0)}ni=1 into a hypothetical
sequence of finite populations under Condition 1.
b−1 (1)(S
bx1,11 − S
bx0,01 ) and γ
b−1 (0)(S
bx0,00 − S
bx1,10 ) are consistent
b1c,RI = S
b0c,RI = S
Corollary 3. γ
xx
xx
b1c,RI − γ
b0c,RI is consistent for βc and follows an asymptotic
for γ1c and γ0c . Furthermore, βbc,RI = γ
normal distribution with covariance matrix


0
S{e0 (0)X} S(εX)
−1 S{e (1)X}
b
cov(βc,RI ) = (πc Sxx,c )
+
−
(πc Sxx,c )−1 ,
n1
n0
n
where we define the residual potential outcomes to be:


 Yi (1) − XiT γ1c ,
 Yi (0) − XiT γ1c ,
0
0
T
Yi (1) − Xi γ0c ,
Yi (0) − XiT γ0c ,
ei (1) =
ei (0) =


T
Yi (1) − Xi γ1c ,
Yi (0) − XiT γ0c ,

Ui = a,
Ui = n,
Ui = c.

(16)

(17)

The idiosyncratic variation is εi = e0i (1) − e0i (0) for unit i, with εi = 0 for Never Takers and
Always Takers, and with εi for Compliers as in (12). The two sets of residuals are not formed from
a regression on all units, but instead the population regression on Compliers alone. As in the ITT
case, we can estimate S{e0 (1)X} and S{e0 (0)X} using their sample analogues; S(εX), however, is
bc1,RI , and for units
unidentifiable. For units with Di = 1, we define the residual eb0i = Yiobs − XiT γ
bc0,RI . Therefore, we can obtain a conservative
with Di = 0, we define the residual eb0i = Yiobs − XiT γ
estimate for the asymptotic covariance (16) by the following sandwich form:
#
"
#
"
0 X)
0 X)
b0 (b
b1 (b
S
e
S
e
−1
−1
−1
b (1) + S
b (0)
b−1 (0).
b (1)
S
S
cov(
c βbc,RI ) = S
xx
xx
xx
xx
n1
n0
As with the ITT analog, so long as we have Assumption 1, randomization itself fully justifies the
theorem and estimators without relying on a model of the observed outcomes.
5.2.2

Two-Stage Least Squares

We now turn to the standard two-stage least squares (TSLS) setting in econometrics (e.g., Angrist
and Pischke, 2008). First, we impose a linear regression model with treatment-covariate interactions:
Yiobs = XiT γ + Di XiT β + ui
17

(i = 1, . . . , n).

Here, the randomness of the observed outcome comes from the randomness of Di and ui . In the
language of econometrics, the treatment received is “endogenous,” i.e., Di and the error term ui
are assumed to be correlated; we therefore use Ti as an instrument for Di . The TSLS estimates
(b
γTSLS , βbTSLS ) are the solutions to the following estimating equations:

n 
X
Xi
−1
bTSLS − Di XiT βbTSLS ) = 0.
(Yiobs − XiT γ
n
Ti Xi

(18)

i=1

This approach is based on M -estimation, though there are many other ways to formalize the
TSLS estimator (e.g., Imbens, 2014). The following theorem shows that the fully-interacted TSLS
estimator βbTSLS is consistent for βc across randomizations.
Theorem 7. Over all randomizations, the TSLS estimator βbTSLS follows an asymptotic normal
distribution with mean βc and covariance matrix


00
S{e00 (0)X} S(εX)
−1 S{e (1)X}
(πc Sxx,c )
+
−
(πc Sxx,c )−1 ,
n1
n0
n
where the residual potential outcomes are defined as


 Yi (0) − XiT (γ∞ + βc ),
 Yi (1) − XiT (γ∞ + βc ),
00
00
T
Yi (1) − Xi γ∞ ,
Yi (0) − XiT γ∞ ,
ei (0) =
ei (1) =


Yi (1) − XiT (γ∞ + βc ),
Yi (0) − XiT γ∞ ,

Ui = a,
Ui = n
Ui = c,

bTSLS , and the idiosyncratic
where γ∞ is the probability limit of the TSLS regression coefficient, γ
treatment effect is εi ≡ e00i (1) − e00i (0).
For variance estimation, define the residual as eb00i = Yiobs − XiT (b
γTSLS + βbTSLS ) for units with
bTSLS for units with Di = 0. We can then use the following sandwich
Di = 1 and eb00i = Yiobs − XiT γ
variance estimator
"
#
"
#
Sb1 (b
e00 X) b−1
Sb0 (b
e00 X) b−1
−1
−1
b
b
b
cov(
c βTSLS ) = Sxx (1)
Sxx (1) + Sxx (0)
Sxx (0),
n1
n0
which has the same probability limit as the Huber–White covariance estimator for βbTSLS . Therefore,
the randomization itself effectively justifies the use of TSLS for estimating systematic treatment
effect variation among Compliers, extending our ITT results.
bTSLS is not, in general, a consistent
Finally, while βbTSLS is a consistent estimator for βc , γ
−1 S
−1
bTSLS converges to γ∞ = Sxx
estimator for γc0 ; that is, γ∞ 6= γc0 . Instead, γ
x0 − πa Sxx Sxx,a βc .
−1 S , the population
In the special case of one-sided noncompliance (i.e., πa = 0), γ∞ = γ0 = Sxx
x0

OLS regression coefficient, among all Compliers and Never Takers, of Y (0) on covariates.
18

5.2.3

Omnibus test for systematic treatment effect variation among Compliers

b for βc , we can use the same Wald-type χ2
With point estimate βb and covariance estimate cov(
c β)
test as in (10) for the presence of systematic treatment effect variation among Compliers. Here, the
estimator can be either randomization-based βbc,RI or TSLS estimator βbTSLS ; the degrees of freedom
are the same, K − 1. Unlike in the ITT case, we are not aware of existing tests for systematic
treatment effect variation among Compliers.

5.3
5.3.1

Idiosyncratic treatment effect variation with noncompliance
Bounding idiosyncratic variation

We now turn to decomposing the overall treatment effect in the presence of noncompliance. In
this setting, we have three sources of treatment effect variation: (i) systematic treatment effect
variation among Compliers, (ii) idiosyncratic treatment effect variation among Compliers, and (iii)
treatment effect variation due to noncompliance.
First, recall that total treatment effect variation is Sτ τ =

Pn

i=1 (τi

− τ )2 /n. We can define a

similar quantity among Compliers:
Sτ τ,c =

n
1 X
I(Ui =c) (τi − τc )2 .
nc
i=1

As in Section 4, we can decompose this variation into systematic and idiosyncratic treatment effect
variation for Compliers, respectively:
Sδδ,c =

n
1 X
I(Ui =c) (δi − τc )2 ,
nc

Sεε,c =

i=1

n
1 X
I(Ui =c) ε2i .
nc
i=1

Because treatment effects for Never Takers and Always Takers are zero, there is no treatment effect
variation for these units. The component of treatment effect variation due to compliance status is
Sτ τ,U =

X

πu (τu − τ )2 .

u=c,a,n

Using τa = τn = 0 and τ = πc τc due to the exclusion restrictions, we have the following theorem
summarizing the relationships among the above components.
Theorem 8. Sτ τ = πc Sτ τ,c + Sτ τ,U , Sτ τ,c = Sδδ,c + Sεε,c , and Sτ τ,U = πc (1 − πc )τc2 .

19

In words, total treatment effect variation has three parts: (i) systematic treatment effect variation among Compliers, πc Sδδ,c ; (ii) idiosyncratic treatment effect variation among Compliers,
πc Sεε,c ; (iii) treatment effect variation due to noncompliance, Sτ τ,U .
As in the ITT case, even though Sεε,c is not identifiable, we can derive bounds in terms of
the marginal distributions of the residuals, {e0i (1) = Yi (1) − XiT γ1c : Ui = c, i = 1, . . . , n} and
{e0i (0) = Yi (0) − XiT γ0c : Ui = c, i = 1, . . . , n}, denoted by F1c (y) and F0c (y), and with marginal
variances, V1c and V0c . Once we estimate these quantities, we can plug them in to Theorem 4
and Corrolary 1 to get our bounds. As compliance status is only partially observed, we have
to estimate these quantities by differencing observed distributions; we defer this and some other
technical details to the Supplementary Material.
5.3.2

Treatment effect decomposition

Since there are two sources of variation—covariates and noncompliance—there are three possible
R2 -type measures. First, we can measure the treatment effect variation explained by noncompliance
alone (i.e., only U ):
2
Rτ,U
=

Sτ τ,U
Sτ τ,U
Sτ τ,U
=
=
.
Sτ τ
Sτ τ,U + πc Sτ τ,c
Sτ τ,U + πc Sδδ,c + πc Sεε,c

Second, we can measure the proportion of treatment effect variation among Compliers explained
by covariates (i.e., only X):
2
Rτ,c
=

Sδδ,c
Sδδ,c
=
.
Sτ τ,c
Sδδ,c + Sεε,c

Third, we can measure the treatment effect variation explained by covariates and noncompliance
(i.e., both X and U ):
2
Rτ,U
X =

Sτ τ,U + πc Sδδ,c
Sτ τ,U + πc Sδδ,c
=
.
Sτ τ
Sτ τ,U + πc Sδδ,c + πc Sεε,c

For each measure, we can use tailored versions of Corollary 1 to construct bounds, or conduct
sensitivity analysis as in Section 4.2, with the sensitivity parameter expressed as the Spearman
correlation between the treatment and control potential outcomes among Compliers.

20

6

Simulation study

6.1

ITT estimators

We simulate completely randomized experiments to evaluate the finite sample performance of the
w , the model-assisted
tests for systematic treatment effect variation based on βbOLS , βbRI , and βbRI

version discussed in the Supplementary Material. Our data generation process is inspired by the
Head Start Impact Study (HSIS) study analyzed in the next section. For a given sample size, we
first generate four independent covariates (X1 , a standard normal, X2 , a binary covariate with
probability 0.5 being 1, X3 , a binary covariate with probability 0.25 being 1, and X4 , a standard
normal). The control potential outcomes are then generated from
Yi (0) = 0.3 + 0.2X1i + 0.3X2i − 0.4X3i + 0.8X4i + ui ,

ui ∼ N (0, σ 2 ).

We select σ 2 = 0.26 to make the marginal variance for the control potential outcomes 1; thus
we can interpret impacts in “effect size” units. The R2 of regressing Y (0) onto the covariates is
approximately 0.74, due to the “pre-test”-like variable X4i . Without X4i , the R2 is about 0.09.
The treatment effects are τi = δi +εi , with (i) either δi = 0.3 for all i, or δi = 0.2+0.1X1i +0.4X3i ;
and (ii) either εi = 0 for all i, or εi ∼ N (0, 0.22 ). All combinations of these two options give the four
cases of (a) no treatment effect variation, (b) only systematic variation, (c) idiosyncratic variation
with no systematic variation, and (d) both systematic and idiosyncratic variation. For an α-level
test of systematic variation, scenarios (a) and (c) should only reject at rate α, while we would like to
see high rejection rates for scenarios (b) and (d). For scenario (d), the Rτ2 is about 0.5; systematic
variation explains a good share of the overall variation.
To generate a synthetic dataset we generated all potential outcomes, randomized units into
treatment with probability 0.6, and then calculated the corresponding observed outcomes. We then
conducted a test for systematic variation using each of our three estimators. For βbRI and βbOLS we
w we also include the fairly predictive X
use X1 , X2 , X3 . For our covariate-adjusted estimator βbRI
4

for adjustment.
Figure 1 shows the power of these tests, with α = 0.05, for different sample sizes. First,
all estimators appear asymptotically valid, consistent with the theoretical results. The OLS and
adjusted estimators are slightly anti-conservative for small n, however, with rejection rates of
around 9%. Second, the OLS estimator appears to have the greatest power in this setting, which is
21

a) none

b) systematic

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

●

●

●

5000

10000

power

●

0.2
●

●

●

●

●

200

500

1000

2000

5000

10000

●

0.0
200

500

1000

2000

0

0.0

0.2
●

●

power

c) idiosyncratic

d) systematic + idiosyncratic

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2
●

0.0
200

●

●

500

1000

N

●

●

●

2000

5000

10000

●

RI

●

●

5000

10000

●

●

●
●

0.0
200

OLS

500

RI (adj)

1000

2000

N

w .
Figure 1: Power of the tests based on βbRI , βbOLS , and βbRI

unsurprising since the true data generating process is a linear model. Finally, covariate adjustment
slightly improves the power of the RI estimator. Overall, in the scenarios we consider, we only
achieve decent levels of power in large samples, although there seems to be reasonable power for
the sample size in the data application, n = 3, 586.

6.2

LATE estimators

We next simulate completely randomized experiments with noncompliance to evaluate the finite
sample performance of the tests for systematic treatment effect variation among Compliers based
on βbc,RI and βbTSLS . We first generated a complete dataset as in the ITT case above, and then
assigned strata membership to all units with probabilities proportional to their covariates. For
Always Takers we then set Yi (0) = Yi (1), and for Never Takers, Yi (1) = Yi (0). The overall ITT is
now reduced to 0.21 (due to the 0 effects of Never Takers and Always Takers), although the CACE

22

is still approximately 0.3. The proportion of Compliers is approximately 68%.
The Compliers have the systematic and idiosyncratic effects described as above. We tested for
the presence of systematic variation for Compliers under the exclusion restrictions. Figure 2 shows
the power of these tests for our RI and TSLS estimators. First, in this scenario, the 2SLS and
the RI estimators are virtually equivalent; the additional adjustment provided by TSLS does not
add significantly to the precision. We see the tests are valid (they even appear conservative) for
cases (a) and (c). Power is reduced compared to the ITT simulation; this is reasonable as power
is effectively a function of the number of Compliers, with additional uncertainty due to partial
information about the identity of Compliers.

power

a) none

b) systematic

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0.0

●

●

●

●

●

●

500

1000

2000

5000

10000

●

●

5000

10000

●

●
●
●

0.0
500

c) idiosyncratic

1000

2000

d) systematic + idiosyncratic

1.0

1.0

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

●

●

5000

10000

power

●

0.0

●

●

●

500

1000

●

2000

N

●

●

5000

10000

●

●
●
●

0.0

RI

500

TSLS

1000

2000

N

Figure 2: Power of the tests based on βbc,RI and βbTSLS .

23

7

Application to the Head Start Impact Study

Established in 1965, Head Start is the largest Federal preschool program in the United States,
serving nearly 1 million low-income three- and four-year-old children each year at a cost of over
$7 billion (Administration for Children and Families, 2015). Researchers and policymakers have
debated Head Start’s effectiveness since its inception, with early randomized trials finding limited
impacts (e.g., Westinghouse Learning Corporation, 1969) and quasi-experimental studies showing
much larger effects (e.g., Currie and Thomas, 1995). Designed in part to settle this debate, the
Head Start Impact Study (HSIS) is a large-scale, nationally representative randomized trial of Head
Start first launched in 2002 (Puma et al., 2010). The Congressional mandate for HSIS included
two broad questions: (1) the program’s overall impact, and (2) how impacts vary across children
and centers. The policy debate has largely focused on this first question; HSIS only found modest
average effects on a range of children’s cognitive and social-emotional outcomes. However, both the
original study and several recent papers argue that these topline results mask important treatment
effect variation (e.g., Bloom and Weiland, 2014; Bitler et al., 2014; Ding et al., 2016; Walters,
2015; Feller et al., 2016). Understanding such variation is critical both for assessing the program’s
benefits and costs and for improving the practice and science of early childhood education.
HSIS collected a rich set of covariates about children and their families, including pre-test score,
child’s age, child’s race, child’s home language, mother’s education level, and mother’s marital
status. At the same time, many potentially important covariates are unavailable. For instance,
while families must be low-income to be eligible for Head Start, HSIS does not include information
on families’ actual income nor other financial details that could be important predictors of program
impact. In addition, Feller et al. (2016) and others argue that that the setting in which a child
would otherwise receive care is an important source of impact variation, although this is not directly
observable.
We now use the methods outlined above to assess treatment effect variation in HSIS. The
original study included n = 4, 400 total children, with n1 = 2, 644 in the treatment group and
n0 = 1, 796 in the control group. Following earlier analyses (Ding et al., 2016) and to simplify
exposition, we restrict our attention to a complete-case subset of the HSIS, with n1 = 2, 238 in the
treatment group and n0 = 1, 348 in the control group (so p1 ≈ 0.62 and p0 ≈ 0.38). Our outcome
of interest is the Peabody Picture Vocabulary Test (PPVT), a widely used measure of cognitive
24

ability in early childhood. To assess treatment effect variation, we consider the full set of childand family-level covariates used in the original HSIS analysis of Puma et al. (2010), including those
mentioned above. After creating dummy variables for factors (e.g., re-coding race), the covariate
matrix has 17 columns. See Figure 3b for a complete list.

7.1

Decomposing variation in the ITT effect

We first explore treatment effect variation for the ITT estimate, beginning with estimating systematic treatment effect variation. We examine three estimators: the randomization-based and OLS
estimators discussed in Section 3, βbRI and βbOLS , and the corresponding model-assisted version of
w . For this latter estimator, we use all
the RI estimator discussed in the Supplementary Material, βbRI

available covariates to adjust the standard estimators, that is, W is the entire vector of covariates.
Omnibus test for systematic treatment effect variation.

We begin by using these estima-

tors for an omnibus test of whether any treatment effect variation is explained by the full set of
w are 0.39 and 0.25,
covariates. The p-values for the unadjusted βbRI estimator and model-assisted βbRI

respectively, which do not show any evidence of treatment effect variation. The OLS estimator,
however, shows much stronger evidence with p = 0.005.
Importantly, all three estimators are based on the same underlying assumptions: the randomization itself justifies all three p-values. And while we expect the unadjusted βbRI to have the lowest
power, it is instructive that the p-value for βbOLS is substantially smaller than the p-value for the
w . As we discuss in Section 3.2, β
bOLS can account for covariate imbalance
covariate-adjusted βbRI

across experimental arms by estimating the Sxx matrix separately for the treatment and control
groups. By contrast, βbRI does not address imbalance in X and instead attempts to residualize out
the Y in order to get a more precise estimate of the relationship of the X to Y for each treatment
arm. Based on the discrepancy in p-values, adjusting for baseline imbalance is clearly important in
this example.
Treatment effect Rτ2 .

Next, we examine how much of the variation could be explained by our

w to estimate the systematic
covariates. Figure 3a shows values of the treatment effect Rτ2 using βbRI

variation. Results are nearly identical using the other estimators. In the worst case of perfect
negative dependence between potential outcomes (not shown), the treatment effect Rτ2 could be as

25

1.0

ITT: Covariates only (X)
LATE: Covariates only (X)
LATE: Covariates and compliance (X + U)

●

0.8

Dual−Language Learner
●
●

●

Child's Race

●

●

0.6

Caregiver's Age

●

Pre−test Score

●

Mother's Marial Status

0.4

Treatment effect R2

●

Mother is Recent Immigrant

●

Mother's Education Level

●

0.2

Child is Male
Both Parents Live with Child

●

0.0

Child's Age ●
●
●
●

0.0

0.2

0.4

0.6

0.8

1.0

0.0

ρ

0.1

0.2

0.3

0.4

0.5

ITT treatment effect R2 with ρ = 1

(a) Overall Rτ2

(b) Rτ2 separately by covariate

Figure 3: Treatment effect Rτ2 , with sensitivity parameter, ρ ∈ [0, 1].
low as 0.01. Assuming that this dependence is nonnegative, the treatment effect Rτ2 ranges from
0.03 to 0.76. While the estimate is clearly sensitive to the unidentifiable sensitivity parameter, the
covariates explain a substantial proportion of treatment effect variation for values of ρ near 1.
We can also use this framework to assess the relative importance of each covariate in terms of
explaining overall treatment effect variation. To do this, we use the model-assisted RI estimator,
w , adjusting for all covariates (i.e., dim(W ) = 17) but restricting systematic treatment effect
βbRI

variation to one covariate at a time. Note that we consider factors (e.g., race) as a group. Figure 3b
shows the resulting estimates for the upper bound of Rτ2 , with lower bound estimates all below 0.01.
Having a mother who is a recent immigrant and dual language learner status (which are highly
correlated in practice) could each explain a substantial proportion of treatment effect variation,
consistent with previous results from Bloom and Weiland (2014) and Bitler et al. (2014). This is
not true for other covariates, like mother’s education level.
Negative correlation between treatment effect and control potential outcomes. Finally, we test whether the individual-level idiosyncratic treatment effects, {εi }ni=1 , are negatively
correlated with the control potential outcomes, {Yi (0)}ni=1 , extending results from Raudenbush and
Bloom (2015). As outlined in the Supplementary Material, we do so by testing whether the variance
26

w : T = 1} is smaller than the variance of {Y obs : T = 0}. This yields a p-value
of {Yiobs − XiT βbRI
i
i
i

of 0.02, which suggests that the unexplained treatment effect is indeed larger for smaller values of
the control potential outcomes. This result is consistent with findings from Bitler et al. (2014) who
use a quantile treatment effect approach.

7.2

Incorporating noncompliance

As with many social experiments, there is substantial noncompliance with random assignment in
HSIS. In the analysis sample we consider here, the estimated proportion of compliance types is
π
bc = 0.69 for Compliers, π
ba = 0.13 for Always Takers, and π
bn = 0.18 for Never Takers. Given the
exclusion restrictions for Always Takers and Never Takers, the treatment effect is therefore zero (by
assumption) for over 30 percent of the sample, suggesting that noncompliance will be an important
component of treatment effect variation.
In the setting with noncompliance, we focus on two estimators for systematic treatment effect
variation among Compliers: the randomization-based estimator, βbc,RI , and the Two-Stage Least
Squares estimator, βbTSLS . We first use these estimators to construct omnibus tests for systematic
treatment effect variation among Compliers. Tests using both estimators show strong evidence for
such variation, with p-value 0.02 using βbc,RI and p-value 0.01 using βbTSLS .
Finally, we turn to decomposing the overall treatment effect. As in the ITT case, we assume that
the potential outcomes have a nonnegative correlation. Figure 3a shows the treatment effect R2
2 = 0.05 to R2 = 0.68. Next, we can calculate treatment
among Compliers, which ranges from Rτ,c
τ,c
2 . In the case of HSIS, this is relatively small—between
effect variation due to noncompliance, Rτ,U

0.01 and 0.16—in part because the overall treatment effect is fairly small. Therefore, the overall
2
treatment effect decomposition due to both covariates and noncompliance, Rτ,U
X , is quite close
2 , as shown in Figure 3a. Taken together, these estimates suggest that there is indeed
to Rτ,c

important treatment effect variation that is neither captured by pre-treatment covariates nor by
noncompliance, consistent with previous results in Ding et al. (2016).

8

Conclusion

In this paper, we propose a broad, flexible framework for assessing and decomposing treatment
effect variation in randomized experiments with and without noncompliance. In general, we believe

27

this is a natural setup for researchers to formulate and investigate a broad range of questions about
impact heterogeneity (e.g., Heckman et al., 1997). Applications include assessing underlying causal
mechanisms and targeting treatments based on individual-level characteristics. Understanding such
variation is also important for the design of experiments. Djebbari and Smith (2008), for example,
argue that characterizing the size of the idiosyncratic treatment effect is useful for determining the
value of additional data collection.
We briefly note several directions for future work. First, our primary purpose was to propose a
framework for analysis rooted in and justified by the randomization itself. As a result, we focused
on the core properties of several relatively simple versions of linear regression and TSLS. We did
not, however, fully explore their practical and finite-sample properties. For example, in future
work, we hope to determine the settings in which model assistance will most improve estimation
and assess the increased power of the OLS approach versus the unbiased RI approach. We are
also investigating how to connect model assisted and OLS approaches to take advantage of both
methods of precision gain. Similarly, there is still much potential improvement in determining ways
of characterizing the degree of heterogeneity, such as with an effect size for the systematic variation.
Second, a natural extension is to use more complex methods to estimate systematic treatment
effects, such as via hierarchical models (Feller and Gelman, 2015) or via machine learning methods (Wager and Athey, 2017), extending the results for the omnibus test and treatment effect
Rτ2 accordingly. While the guarantees from randomization are clearly weaker in such settings,
researchers can assess these tradeoffs themselves. For example, hierarchical modeling would be
especially useful in the Head Start Impact Study due to the multi-site design (Bloom and Weiland,
2014).
Third, a question of increasing practical importance is the generalizability of experimental
results to a given target population (Stuart et al., 2011). We believe that the treatment effect
Rτ2 is a critical measure for assessing the credibility of these generalizations. In short, if there is
substantial idiosyncratic treatment effect variation, i.e., Rτ2 is small, then researchers should be
wary of using observed covariates to extrapolate treatment effects.
Finally, a question is how to extend this treatment effect variation framework to non-randomized
settings. While the results would necessarily rest on much stronger assumptions, many settings
already use an as-if-randomized framework, such as in observational studies (Rosenbaum, 2002;

28

Imbens and Rubin, 2015). Under this approach, extensions should be natural.

References
A. Abadie. Semiparametric instrumental variable estimation of treatment response models. Journal
of Econometrics, 113:231–263, 2003.
Administration for Children and Families. Head Start program facts, fiscal year 2014. Available at
https://eclkc.ohs.acf.hhs.gov/hslc/data/factsheets/docs/hs-program-fact-sheet-2014.pdf, 2015.
J. D. Angrist and J. Pischke. Mostly Harmless Econometrics: An Empiricist’s Companion. Princeton: Princeton University Press, 2008.
J. D. Angrist, G. W. Imbens, and D. B. Rubin. Identification of causal effects using instrumental
variables. Journal of the American Statistical Association, 91:444–455, 1996.
J. D. Angrist, P. A. Pathak, and C. R. Walters. Explaining charter school effectiveness. American
Economic Journal: Applied Economics, 5(4):1–27, 2013.
P. M. Aronow, D. P. Green, and D. K. Lee. Sharp bounds on the variance in randomized experiments. The Annals of Statistics, 42:850–871, 2014.
S. Athey and G. Imbens. Recursive partitioning for heterogeneous causal effects. Proceedings of
the National Academy of Sciences, 113(27):7353–7360, 2016.
A. Berrington de González and D. R. Cox. Interpretation of interaction: A review. The Annals of
Applied Statistics, 1:371–385, 2007.
M. Bitler, H. Hoynes, and T. Domina. Experimental Evidence on Distributional Effects of Head
Start. Working Paper, 2014.
A. S. Blinder. Wage discrimination: reduced form and structural estimates. Journal of Human
resources, 8:436–455, 1973.
H. S. Bloom and C. Weiland. To what extent do the effects of Head Start on enrolled children vary
across sites? Working Paper, 2014.

29

J. Bound, D. A. Jaeger, and R. M. Baker. Problems with instrumental variables estimation when the
correlation between the instruments and the endogenous explanatory variable is weak. Journal
of the American statistical association, 90:443–450, 1995.
W. G. Cochran. Sampling Techniques. New York: John Wiley & Sons, 3rd edition, 1977.
D. R. Cox. Interaction (with discussion). International Statistical Review, 52:1–24, 1984.
R. K. Crump, V. J. Hotz, G. W. Imbens, and O. A. Mitnik. Nonparametric tests for treatment
effect heterogeneity. Review of Economics and Statistics, 90:389–405, 2008.
J. Currie and D. Thomas. Does Head Start make a difference? American Economic Review, 85(3):
341–364, 1995.
P. Ding. A paradox from randomization-based causal inference. arXiv preprint arXiv:1402.0142,
2014.
P. Ding, A. Feller, and L. W. Miratrix. Randomization inference for treatment effect variation.
Journal of the Royal Statistical Society, Series B (Statistical Methodology), 78:655–671, 2016.
H. Djebbari and J. Smith. Heterogeneous impacts in PROGRESA. Journal of Econometrics, 145:
64–80, 2008.
A. Feller and A. Gelman. Hierarchical models for causal effects. Emerging Trends in the Social and
Behavioral Sciences: An Interdisciplinary, Searchable, and Linkable Resource, 2015.
A. Feller, T. Grindal, L. Miratrix, and L. C. Page. Compared to what? variation in the impacts
of early childhood education by alternative care type. The Annals of Applied Statistics, 10(3):
1245–1285, 2016.
R. A. Fisher. The Design of Experiments. Edinburgh: Oliver & Boyd, 1st edition, 1935.
C. B. Fogarty. Regression assisted inference for the average treatment effect in paired experiments.
arXiv preprint arXiv:1612.05179, 2016.
C. E. Frangakis and D. B. Rubin. Principal stratification in causal inference. Biometrics, 58:21–29,
2002.

30

M. Fréchet. Sur les tableaux de corrélation dont les marges son données. Annals Universite de
Lyon, Sect. A. Ser. 3, 14:53–77, 1951.
D. P. Green and H. L. Kern. Modeling Heterogeneous Treatment Effects in Survey Experiments
with Bayesian Additive Regression Trees. The Public Opinion Quarterly, 76:491–511, 2012.
J. Hájek. Limiting distributions in simple random sampling from a finite population. Publications
of the Mathematics Institute of the Hungarian Academy of Science, 5:361–74, 1960.
J. J. Heckman, J. Smith, and N. Clements. Making the most out of programme evaluations
and social experiments: Accounting for heterogeneity in programme impacts. The Review of
Economic Studies, 64:487–535, 1997.
J. L. Hill. Bayesian nonparametric modeling for causal inference. Journal of Computational and
Graphical Statistics, 20:217–240, 2011.
W. Hoeffding. Masstabinvariante korrelationsmasse für diskontinuierliche verteilungen. Arkiv fr
matematischen Wirtschaften und Sozialforschung, 7:49–70, 1941.
Y. Huang, P. B. Gilbert, and H. Janes. Assessing treatment-selection markers using a potential
outcomes framework. Biometrics, 68:687–696, 2012.
P. J. Huber. The behavior of maximum likelihood estimates under nonstandard conditions. In Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, volume 1,
pages 221–233, 1967.
K. Imai and M. Ratkovic. Estimating treatment effect heterogeneity in randomized program evaluation. The Annals of Applied Statistics, 7:443–470, 2013.
G. Imbens. Instrumental variables: An econometrician’s perspective (with discussion). Statistical
Science, 29:323–358, 2014.
G. W. Imbens and D. B. Rubin. Causal Inference in Statistics, and in the Social and Biomedical
Sciences. New York: Cambridge University Press, 2015.
O. Kempthorne. The Design and Analysis of Experiments. New York: Wiley, 1952.
E. L. Lehmann. Elements of Large-Sample Theory. New York: Springer, 1998.
31

X. Li and P. Ding. General forms of finite population central limit theorems with applications to
causal inference. Journal of the American Statistical Association, page in press, 2017.
W. Lin. Agnostic notes on regression adjustments to experimental data: reexamining Freedman’s
critique. The Annals of Applied Statistics, 7:295–318, 2013.
R. A. Matsouaka, J. Li, and T. Cai. Evaluating marker-guided treatment selection strategies.
Biometrics, 70:489–499, 2014.
J. A. Middleton and P. M. Aronow. Unbiased estimation of the average treatment effect in clusterrandomized experiments. Statistics, Politics and Policy, 6:39–75, 2015.
R. B. Nelsen. An Introduction to Copulas. New York: Springer, 2nd edition, 2007.
J. Neyman. On the application of probability theory to agricultural experiments. Essay on principles. Section 9. Statistical Science, 5:465–472, 1923.
R. Oaxaca. Male-female wage differentials in urban labor markets. International Economic Review,
14:693–709, 1973.
M. Puma, S. Bell, R. Cook, C. Heid, G. Shapiro, P. Broene, F. Jenkins, P. Fletcher, L. Quinn,
J. Friedman, et al. Head start impact study: Final report. Technical report, Department of
Health and Human Services, Administration for Children and Families, Washington DC, 2010.
S. W. Raudenbush and H. S. Bloom. Learning about and from a distribution of program impacts
using multisite trials. American Journal of Evaluation, 36(4):475–499, 2015.
P. R. Rosenbaum. Reduced sensitivity to hidden bias at upper quantiles in observational studies
with dilated treatment effects. Biometrics, 55:560–564, 1999.
P. R. Rosenbaum. Observational Studies. New York: Springer, 2nd edition, 2002.
P. R. Rosenbaum. Confidence intervals for uncommon but dramatic responses to treatment. Biometrics, 63:1164–1171, 2007.
D. B. Rubin. Estimating causal effects of treatments in randomized and nonrandomized studies.
Journal of Educational Psychology, 66:688–701, 1974.

32

D. B. Rubin. Comment on “Randomization analysis of experimental data: the Fisher randomization
test” by D. Basu. Journal of the American Statistical Association, 75:591–593, 1980.
C.-E. Särndal, B. Swensson, and J. Wretman.

Model-Assisted Survey Sampling.

New York:

Springer, 2003.
D. O. Staiger and J. H. Stock. Instrumental variables regression with weak instruments. Econometrica, 65:557–586, 1997.
E. A. Stuart, S. R. Cole, C. P. Bradshaw, and P. J. Leaf. The use of propensity scores to assess
the generalizability of results from randomized trials. Journal of the Royal Statistical Society:
Series A (Statistics in Society), 174(2):369–386, 2011.
S. Wager and S. Athey. Estimation and inference of heterogeneous treatment effects using random
forests. Journal of the American Statistical Association, (just-accepted), 2017.
C. R. Walters. Inputs in the production of early childhood human capital: Evidence from head
start. American Economic Journal: Applied Economics, 7(4):76–102, 2015.
Westinghouse Learning Corporation. The Impact of Head Start: An Evaluation of the Effects of
Head Start on Children’s Cognitive and Affective Development, Volume 1: Report to the Office of
Economic Opportunity. Athens, Ohio: Westinghouse Learning Corporation and Ohio University,
1969.
H. White. A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity. Econometrica, pages 817–838, 1980.
A. Zhao, P. Ding, R. Mukerjee, and T. Dasgupta. Randomization-based causal inference from
split-plot designs. Annals of Statistics, page in press, 2017.

33

Supplementary Material
for
“Decomposing Treatment Effect Variation”
Appendix A gives all the proofs and Appendix B provides the additional commentary mentioned
in the main text. The finite population central limit theorem (FPCLT) we use for our asymptotic
proofs is Theorem 5 of Li and Ding (2017), which requires some mild moment conditions on the
covariates and potential outcomes, as outlined in the main text.

Appendix A

Lemmas and Proofs

Before we prove Theorem 1, we provide a few lemmas to ease the notational burden and amount of
algebra of subsequent calculations. These lemmas allow us to derive expressions for our estimators
in terms of matrix algebra rather than the summation-style approach typically seen for Neymanstyle derivations in the literature.
To begin, let 1n = (1, . . . , 1)T and 0n = (0, . . . , 0)T be column vectors of length n, and In be the
n × n identity matrix. Then Sn = In − n−1 1n 1Tn is the projection matrix orthogonal to 1n with
Sn 1n = 0n . Under this formulation, the covariance matrix of the treatment assignment vector is a
scaled projection matrix orthogonal to 1n , as shown in the following lemma.
Lemma A.1. The treatment assignment vector T of a completely randomized experiment has
E(T ) =

n1
1n ,
n

cov(T ) =

n1 n0
Sn .
n(n − 1)

Proof of Lemma A.1. The conclusions follow from
E(Ti ) =

n1
,
n

var(Ti ) =

n1 n0
,
n2

cov(Ti , Tj ) = −

n1 n0
,
2
n (n − 1)

(i 6= j).

The projection matrix Sn acts as a covariance operator as illustrated by the following lemma.
Lemma A.2. Let Ui , Vi ∈ RK be column vectors of length K. Define U = [U1 , U2 , . . . , Un ] and
P
V = [V1 , V2 , . . . , Vn ] ∈ RK×n as two matrices of dimension K × n. If Ū = n−1 ni=1 Ui = n−1 U1
P
and V̄ = n−1 ni=1 Vi = n−1 V1, then
USn V =
T

n
X

(Ui − Ū )(Vi − V̄ )T .

i=1

1

(A.1)

In particular, when Ui = Vi ,
VSn V T =

n
X
(Vi − V̄ )(Vi − V̄ )T = (n − 1)S(V ).
i=1

Proof of Lemma A.2. The left hand side of (A.1) is equal to
USn V T = UV T − n−1 (U1n ) (V1n )T =

n
X

Ui ViT − n−1 (nŪ )(nV̄ )T =

i=1

n
X

Ui ViT − nŪ V̄ T ,

i=1

which is the same as the right hand side of (A.1).
Theorem 1: A generalized, vector-outcome version of Neyman. To prove the generalized
Neyman result, we bundle our vector potential outcomes into matrices and use the above lemmas
to obtain their covariance matrix. The theorem is exact, no asymptotics. Using the FPCLT to
show that the estimator has an approximately Normal distribution, allowing for classic testing and
inference, is a separate, subsequent step.
Proof of Theorem 1. Define V1 = [V1 (1), . . . , Vn (1)] and V0 = [V1 (0), . . . , Vn (0)] as the matrices of
the potential outcomes. Then the Neymanian simple difference in means estimator has the following
representation:
τbV

= V̄1obs − V̄0obs
n
n
1 X
1 X
=
Ti Vi (1) −
(1 − Ti )Vi (0)
n1
n0
i=1

i=1

1
1
V1 T − V0 (1 − T )
=
n
n0
1

1
V1 V0
=
+
T − V0 1.
n1 n0
n0
Now the unbiasedness of τbV follows from the linearity of the expectation and Lemma A.1. For
the covariance, note the second term in the above is constant, and so is not involved. Applying
Lemmas A.1 and A.2, we can obtain the covariance matrix of τbV :




V1 V0 T
V1 V0
cov(b
τV ) =
+
cov(T )
+
n1 n0
n1 n0




n1 n0
V1 V0
V1 V0 T
=
+
Sn
+
n(n − 1) n1 n0
n1 n0


n1 n0
1
1
1
1
T
T
T
T
=
V1 Sn V1 + 2 V0 Sn V0 +
V0 Sn V1 +
V1 Sn V0
n(n − 1) n21
n1 n0
n1 n0
n0
n0
n1
1
=
S{V (1)} +
S{V (0)} +
(V0 Sn V1T + V1 Sn V0T ).
nn1
nn0
n(n − 1)
2

To simplify the third term, we use the fact abT + baT = aaT + bbT − (a − b)(a − b)T for two column
vectors a and b, we have
{Vi (1) − V̄ (1)}{Vi (0) − V̄ (0)}T + {Vi (0) − V̄ (0)}{Vi (1) − V̄ (1)}T
= {Vi (1) − V̄ (1)}{Vi (1) − V̄ (1)}T + {Vi (1) − V̄ (1)}{Vi (1) − V̄ (1)}T
−{Vi (1) − Vi (0) − V̄ (1) + V̄ (0)}{Vi (1) − Vi (0) − V̄ (1) + V̄ (0)}T .
Summing over i = 1, . . . , n and applying Lemma A.2, we have
V0 Sn V1T
V1 Sn V0T
+
= S{V (1)} + S{V (0)} − S{V (1) − V (0)}.
n−1
n−1
Therefore, the covariance of τbV can be simplified as:
cov(b
τV ) =
=

n1
1
n0
S{V (1)} +
S{V (0)} + [S{V (1)} + S{V (0)} − S{V (1) − V (0)}]
nn1
nn0
n
S{V (1)} S{V (0)} S{V (1) − V (0)}
+
−
.
n1
n0
n

Theorem 2: Behavior of βbRI .

To show properties of βbRI we express the systematic variation

as a vector of new potential outcomes of the original outcome scaled by the different covariates of
interest. This allows for immediate use of Theorem 1.
bxt is the sample mean for {Xi Y obs : Ti = t, i = 1, . . . , n} =
Proof of Theorem 2. Because S
i
{Xi Yi (t) : Ti = t, i = 1, . . . , n}, it is unbiased for the population mean Sxt . Thus, the estima−1 is fixed and the expectation is linear. Its sampling covariance
tor βbRI is also unbiased for β as Sxx

over all possible randomizations is
−1
bx1 − S
bx0 )S −1 .
cov(βbRI ) = Sxx
cov(S
xx

Therefore, we need only to obtain the covariance of
n
n
1 X
1 X
obs
b
b
Sx1 − Sx0 =
Ti Xi Yi −
(1 − Ti )Xi Yiobs ,
n1
n0
i=1

i=1

which is the difference between the sample means of {Xi Yi (1) : i = 1, . . . , n} and {Xi Yi (0) : i = 1, . . . , N }
under treatment and control. Viewing Xi Yiobs as a vector outcome in a completely randomized
experiment, we can apply Theorem 1 to obtain
bx1 − S
bx0 ) = S{XY (1)} + S{XY (0)} − S(Xτ ) ,
cov(S
n1
n0
n
3

which completes the proof.
Theorem 3: Behavior of βbOLS .

We first use the well-known fact that the estimate from a

OLS model with treatment fully interacted with covariates is equivalent to separate regressions of
bOLS
outcome onto covariates for the control and treatment groups. This means we can obtain γ
by running a regression of Y obs onto X using the control group data, and γ\
+ β OLS by running
regression of Y obs onto X using the treatment group data, giving estimated coefficients of
b−1 S
b
bOLS = S
γ
xx,0 x0
and
b−1 S
b
b−1 b
βbOLS = S
xx,1 x1 − Sxx,0 Sx0 .
As a quick heuristic argument for this, consider that the maximization problem for the interacted
model will separate into two components, one for each group. Then re-parameterize to get the
above.
We now prove the properties of βbOLS . Here we have to use asymptotics for the entire theorem,
unlike the case of βbRI , where the mean and covariance are exact and the asymptotics are only
needed for the asymptotic normality of the estimator.
Proof of Theorem 3. First expand the difference of βbOLS and β as
bx0 − S
bxx,0 γ0 ),
bx1 − S
bxx,1 γ1 ) − S
b−1 (S
b−1 (S
βbOLS − β = S
xx,0
xx,1
This will be close to the related quantity of
−1 b
bxx,1 γ1 ) − S −1 (S
bx0 − S
bxx,0 γ0 ).
∆ = Sxx
(Sx1 − S
xx

(A.2)

For the above to make sense and hold, we here need our asymptotic framework. In particular,
we need the associated moment conditions described in the main text. We next observe that the
difference between βbOLS − β and ∆ is of higher order, because
b−1 − S −1 )(S
bx1 − S
bxx,1 γ1 ) − (S
b−1 − S −1 )(S
bx0 − S
bxx,0 γ0 ) (A.3)
(βbOLS − β) − ∆ = (S
xx
xx
xx,1
xx,0
= OP (n−1/2 )OP (n−1/2 ) − OP (n−1/2 )OP (n−1/2 ) = OP (n−1 ),

(A.4)

following from the FPCLT for the four terms in (A.3). This is an argument commonly used in the
survey sampling literature for ratio estimators (Cochran, 1977).
4

We next focus on the asymptotic distribution of ∆, because the asymptotic distribution of
βbOLS − β will be the same. Further simplify (A.2) as
#
"
n
n
X
X
1
1
−1
Ti Xi ei (1) −
(1 − Ti )Xi ei (0) ,
∆ = Sxx
n1
n0
i=1

(A.5)

i=1

where ei (1) = Yi (1) − XiT γ1 and ei (0) = Yi (0) − XiT γ0 are the residual potential outcomes. (To see
bx1 and S
bxx,1 are sums over the treatment units, and we
the above, note, for example, that both S
can factor out an Xi to get Xi times the difference in the Yi and predicted Yi .)
Applying Theorem 1 to the vector outcome Xe, we obtain the covariance matrix of ∆, to which
βbOLS − β converges to due to (A.4). The asymptotic normality follows from the representation
(A.5) and the FPCLT.
Theorem 4: Bounds for Rτ2 .

To prove Theorem 4, we need to invoke the following Fréchet–

Hoeffding inequality (Hoeffding, 1941; Fréchet, 1951; Heckman et al., 1997; Aronow et al., 2014).
Lemma A.3. If we know only the marginal distributions of two random variables X ∼ FX (x) and
Y ∼ FY (y), then E(XY ) can be sharply bounded by
Z
Z 1
−1
−1
FX (u)FY (1 − u)du ≤ E(XY ) ≤

0

0

1

FX−1 (u)FY−1 (u)du.

Lemma A.3 immediately implies the following bound for var(X − Y ) if E(X − Y ) = 0.
Lemma A.4. If we know only the marginal distributions X ∼ FX (x), Y ∼ FY (y) and E(X − Y ) =
0, then var(X − Y ) can be sharply bounded by
Z 1
Z 1
−1
−1
2
{FX−1 (u) − FY−1 (1 − u)}2 du
{FX (u) − FY (u)} du ≤ var(X − Y ) ≤
0

0

Proof of Lemma A.4. The variance var(X − Y ) can be decomposed as
var(X − Y ) = E(X − Y )2 = E(X 2 ) + E(Y 2 ) − 2E(XY ),
which depends on the following three terms:
2

E(X ) =

R

x2 dFX (x)

Z

1

=
0

2

E(Y ) =
Z
0

1

FX−1 (u)FY−1 (1 − u)du ≤

R1

−1
2
0 {FY (u)} du

Z
0

Z
E(XY )

1

=
≤
0

1

{FX−1 (u)]}2 du,
{FY−1 (1 − u)}2 du,
FX−1 (u)FY−1 (u)du.

Plug the above expressions into the variance of X − Y to obtain the desired bounds.
5

Applying Lemma A.4, we can easily prove Theorem 4.
Proof of Theorem 4. Because Sτ τ = Sδδ + Sεε , we need only to bound Sεε , which is the finite
population variance of εi = {Yi (1) − XiT γ1 } − {Yi (0) − XiT γ0 } = ei (1) − ei (0). We can identify
the marginal distributions of {ei (1) : i = 1, . . . , n} and {ei (0) : i = 1, . . . , n}, and also know that
P
n−1 ni=1 εi = 0. Therefore, the bounds in Lemma A.4 imply the bounds in Theorem 4.
Theorem 5: Sensitivity analysis.
Proof of Theorem 5. The joint distribution of (U1 , U0 ) is
C(u1 , u0 ) = P (U1 ≤ u1 , U0 ≤ u0 )
= ρP (U0 ≤ u1 , U0 ≤ u0 ) + (1 − ρ)P (V0 ≤ u1 , U0 ≤ u0 )
= ρ min(u1 , u0 ) + (1 − ρ)u1 u0 .
Therefore, the distribution function C(u1 , u0 ) is a weighted average of min(u1 , u0 ) = CR (u1 , u0 )
and u1 u0 = CI (u1 , u0 ), i.e., the joint distributions when U1 = U0 and U1 U0 , respectively.
According to Nelsen (2007, Theorem 5.1.6), Spearman’s rank correlation coefficient between
e(1) and e(0) is
Z

1Z 1

1Z 1

Z
{C(u1 , u0 ) − u1 u0 }du1 du0 = 12ρ

12
0

{min(u1 , u0 ) − u1 u0 }du1 du0
0 Z 01

Z u1
1
= 12ρ 2
du1
u0 du0 −
4
0
0
= 12ρ(1/3 − 1/4) = ρ.

0

To complete the proof of the theorem, we need only to show that the covariance between e(1)
and e(0) is linear in ρ, which follows from
1Z 1

Z
0

0

Z

F1−1 (u1 )F0−1 (u0 )dC(u1 , u0 )

1Z 1

= ρ
0

Z
= ρ
0

0
1

F1−1 (u1 )F0−1 (u0 )dCR (u1 , u0 )

F1−1 (u)F0−1 (u)du + (1 − ρ)

Z

1

0

6

Z

1Z 1

+ρ
0

F1−1 (u)du

0

Z
0

F1−1 (u1 )F0−1 (u0 )dCI (u1 , u0 )
1

F0−1 (u)du.

Theorem 6: Extending to non-compliance.

Theorem 6 shows how to estimate the outcome-

to-covariate relationships of the Compliers by estimating different aggregate covariance relationships
across all the strata for different observed groups and then taking differences. Due to the exclusion
restriction for the Never Takers and Always Takers, this gives our desired relationships for the
Compliers only.
First, a small bit of notation of, due to the exclusion restrictions for Never Takers and Always
Takers, defining the population covariance between X and Y (1) = Y (0) within stratum U = a and
U = n as
Sx.,u

n
n
1 X
1 X
=
I(Ui =u) Xi Yi (1) =
I(Ui =u) Xi YI (0),
nu
nu
i=1

(u = a, n).

i=1

Proof of Theorem 6. We first create an estimator for Sxx,c . From the observed data with (Ti , Di ) =
(1, 1), we have
)
)
(
(
n
n
n
1 X
1 X
1 X
T
T
T
Ti Di Xi Xi
Ti I(Ui =a) Xi Xi +
Ti I(Ui =c) Xi Xi
= E
E
n1
n1
n1
i=1

i=1

i=1

= πa Sxx,a + πc Sxx,c .

(A.6)

Similar to (A.6), we have
)
n
1 X
T
Ti (1 − Di )Xi Xi
= πn Sxx,n ,
E
n1
i=1
(
)
n
1 X
E
(1 − Ti )Di Xi XiT
= πa Sxx,a ,
n0
i=1
(
)
n
1 X
E
(1 − Ti )(1 − Di )Xi XiT
= πn Sxx,n + πc Sxx,c .
n0
(

(A.7)

(A.8)

(A.9)

i=1

Subtracting the left sides of (A.8) from (A.6), or subtracting the left sides of (A.7) from (A.9), give
unbiased estimators for πc Sxx,c .
Second, analogous to the Sxx,c , we consider the sample covariances between X and Y obs to
obtain estimators for Sx1,c and Sx0,c . From the observed data with (Ti , Di ) = (1, 1), we have
)
(
)
(
n
n
n
1 X
1 X
1 X
obs
Ti Di Xi Yi
= E
Ti I(Ui =a) Xi Yi (1) +
Ti I(Ui =c) Xi Yi (1)
E
n1
n1
n1
i=1

i=1

i=1

= πa Sx.,a + πc Sx1,c .

7

(A.10)

Similar to (A.10), we have
)
n
1 X
Ti (1 − Di )Xi Yiobs = πn Sx.,n ,
E
n1
i=1
(
)
n
1 X
E
(1 − Ti )Di Xi Yiobs = πa Sx.,a ,
n0
i=1
)
(

(
E

n
1 X
(1 − Ti )(1 − Di )Xi Yiobs
n0

= πn Sx.,n + πc Sx0,c .

(A.11)

(A.12)

(A.13)

i=1

Subtracting (A.12) from (A.10), and subtracting (A.11) from (A.13), we obtain the results in
(15).
Corollary 3: Behavior of βbc,RI .

Theorem 6 shows how to obtain unbiased estimates of the

components of our estimator, which we can then plug in to obtain a consistent estimator of βc . We
next show how this plug-in estimator behaves.
Proof of Corollary 3. First we write
bxx,11 − S
bxx,01 )−1 {S
bx1,11 − S
bx0,01 − (S
bxx,11 − S
bxx,01 )γ1c }
βbc,RI − βc = (S
bxx,00 − S
bxx,10 )−1 {S
bx0,00 − S
bx1,10 − (S
bxx,00 − S
bxx,10 )γ0c },
−(S
second we introduce
bx1,11 − S
bx0,01 − (S
bxx,11 − S
bxx,01 )γ1c }
∆c = (πc Sxx,c )−1 {S
bx0,00 − S
bx1,10 − (S
bxx,00 − S
bxx,10 )γ0c },
−(πc Sxx,c )−1 {S
third we observed that the difference between βbc,RI − βc and ∆c has higher order following the
same argument as (A.4). Therefore, we need only to find the asymptotic distribution of ∆c .

8

Simple algebra gives
∆c = (πc Sxx,c )−1

n
n
h1 X
1 X
Ti Di Xi Yi (1) −
(1 − Ti )Di Xi Yi (0)
n1
n0
i=1

−
−
+

1
n1
1
n0
1
n0

n
X
i=1
n
X
i=1
n
X

Ti Di Xi XiT γc1 +

i=1

1
n0

n
X

(1 − Ti )Di Xi XiT γc1

i=1

(1 − Ti )(1 − Di )Xi Yi (0) +

n
1 X
Ti (1 − Di )Xi Yi (1)
n1

(1 − Ti )(1 − Di )Xi XiT γc0 −

i=1

i=1
n
X

1
n1

Ti (1 − Di )Xi XiT γc0

i

i=1

n
n
n
h1 X
1 X
1 X
= (πc Sxx,c )−1
Ti I(Ui =a) Xi Yi (1) +
Ti I(Ui =c) Xi Yi (1) −
(1 − Ti )I(Ui =a) Xi Yi (0)
n1
n1
n0
i=1

i=1

i=1

n
n
n
1 X
1 X
1 X
−
Ti I(Ui =a) Xi XiT γc1 −
Ti I(Ui =c) Xi XiT γc1 +
(1 − Ti )I(Ui =a) Xi XiT γc1
n1
n1
n0
i=1

i=1

i=1

n
n
n
1 X
1 X
1 X
−
(1 − Ti )I(Ui =n) Xi Yi (0) −
(1 − Ti )I(Ui =c) Xi Yi (0) +
Ti I(Ui =n) Xi Yi (1)
n0
n0
n1
i=1

i=1

i=1

n
n
n
i
1 X
1 X
1 X
T
T
+
(1 − Ti )I(Ui =n) Xi Xi γc0 +
(1 − Ti )I(Ui =c) Xi Xi γc0 −
Ti I(Ui =n) Xi XiT γc0
n0
n0
n1
i=1

i=1

i=1

n
n X


−1 1
Ti Xi I(Ui =a) (Yi (1) − XiT γc1 ) + I(Ui =n) (Yi (1) − XiT γc0 ) + I(Ui =c) (Yi (1) − XiT γc1 )
= (πc Sxx,c )
n1
i=1

n

o
1 X
(1 − Ti )Xi I(Ui =a) (Yi (0) − XiT γc1 ) + I(Ui =n) (Yi (0) − XiT γc0 ) + I(Ui =c) (Yi (0) − XiT γc0 ) .
−
n0
i=1

According to the definitions of the residual potential outcomes e0i (1) and e0i (0) in the main text, the
above formula reduces to
"
βec,RI − βc = (πc Sxx,c )

−1

#
n
n
1 X
1 X
0
0
Ti Xi ei (1) −
(1 − Ti )Xi ei (0) .
n1
n0
i=1

(A.14)

i=1

The representation in (A.14) implies the asymptotic covariance matrix according to Theorem 1 and
the asymptotic normality of βec,RI according to the FPCLT.
Theorem 7: Behavior of βbTSLS .

While the amount of notation and matrix algebra is consider-

ably more in scope, the overall structure of the proof follows the earlier one for the OLS estimator
for the ITT. In particular, we show the estimator asymptotically converges to a more tractable
version that has a fixed portion, and then use the usual covariance argument on the remaining

9

terms. Before doing this, we first show the probability limits of the estimator by working through
the matrix algebra.
Proof of Theorem 7. First, we find the probability limits of the TSLS estimators:
)−1 ( n 
)
( n 




bTSLS
γ
1 X Xi
1 X Xi
T
T
obs
(Xi , Di Xi )
Yi
=
Ti Xi
Ti Xi
βbTSLS
n
n
i=1
i=1
 −1 Pn
−1  −1 Pn

P
n P i=1 Xi XiT
n−1P ni=1 Di Xi XiT
n P i=1 Xi Yiobs
=
n−1 ni=1 Ti Xi XiT n−1 ni=1 Ti Di Xi XiT
n−1 ni=1 Ti Xi Yiobs

−1  
A B
G
P
−→
.
(A.15)
C D
H
The above term A is A = Sxx , and terms (B, C, D, G, H) are the population limits of the sample
quantities. We will find each of them. Term B is
)
)
( n
( n
n
X
X
1
1
1X
Di Xi XiT = E
Ti Di Xi XiT +
(1 − Ti )Di Xi XiT
B = E
n
n
n
i=1
i=1
i=1
( n
)
n
n
X
X
X
1
1
1
Ti I(Ui =a) Xi XiT +
Ti I(Ui =c) Xi XiT +
(1 − Ti )I(Ui =a) Xi XiT
= E
n
n
n
i=1

i=1

i=1

= p1 πa Sxx,a + p1 πc Sxx,c + p0 πa Sxx,a
= πa Sxx,a + p1 πc Sxx,c .

P
Term C is C = E n−1 ni=1 Ti Xi XiT = p1 Sxx . Term D is
( n
)
( n
)
n
1X
1X
1X
T
T
T
D = E
Ti Di Xi Xi = E
Ti I(Ui =a) Xi Xi +
Ti I(Ui =c) Xi Xi
n
n
n
i=1

i=1

i=1

= p1 πa Sxx,a + p1 πc Sxx,c .
Term G is
(
G=E

n

1X
Xi Yiobs
n

)

(
=E

i=1

n

n

i=1

i=1

1X
1X
Ti Xi Yiobs +
(1 − Ti )Xi Yiobs
n
n

)
= p1 Sx1 + p0 Sx0 .


P
Term H is H = E n−1 ni=1 Ti Xi Yiobs = p1 Sx1 . We apply the following formula for the inverse
of a block matrix:


A B
C D

−1


−1
−1
SD
−A−1 BSA
,
−1
−1
−D −1 CSD
SA


=

where SD = A − BD −1 C and SA = D − CA−1 B are the Schur complements of blocks D and
A. Omitting some tedious matrix algebra, we obtain
SD = p0 πc Sxx,c (πa Sxx,a + πc Sxx,c )−1 Sxx ,
10

SA = p1 p0 πc Sxx,c ,

and the inverse of the block matrix is

−1  −1 −1 −1

−1
−1 (π S
−1
p0 πc Sxx (πa Sxx,a + πc Sxx,c )Sxx,c
−p−1
p−1
πc−1 Sxx
A B
a xx,a + p1 πc Sxx,c )Sxx,c
1
0
=
.
−1 −1 −1
−1 −1
C D
−p−1
p−1
0 πc Sxx,c
1 p0 πc Sxx,c
bTSLS is
Therefore, according to (A.15), the probability limit of γ
−1 −1 −1 −1
−1 −1
−1
−1
p−1
0 πc Sxx (πa Sxx,a + πc Sxx,c )Sxx,c (p1 Sx1 + p0 Sx0 ) − p1 p0 πc Sxx (πa Sxx,a + p1 πc Sxx,c )Sxx,c (p1 Sx1 )
−1
−1
−1
= Sxx
Sx0 − πa πc−1 Sxx
Sxx,a Sxx,c
(Sx1 − Sx0 )
−1
= γ0 − πa Sxx
Sxx,a βc ≡ γ∞ ,

(A.16)

and the probability limit of βbTSLS is
−1 −1 −1 −1
−1 −1
−1 −1
−p−1
0 πc Sxx,c (p1 Sx1 + p0 Sx0 ) + p1 p0 πc Sxx,c (p1 Sx1 ) = πc Sxx,c (Sx1 − Sx0 ) = βc , (A.17)

where we use Sx1 − Sx0 = πc (Sx1,c − Sx0,c ), which is guaranteed by exclusion restrictions.
We next find the asymptotic distribution of βbTSLS . Following the derivation in Corollary 3, we
first write
)−1 ( n 

   ( X
)

n 
bTSLS
γ
1
1 X Xi (Yiobs − XiT γ∞ − Di XiT βc )
γ∞
Xi
T
T
−
=
,
(Xi , Di Xi )
βc
Ti Xi
Ti Xi (Yiobs − XiT γ∞ − Di XiT βc )
βbTSLS
n
n
i=1

i=1

then introduce


∆TSLS

A
=
C

A
=
C

A
=
C

)
−1 ( X
n 
1
Xi (Yiobs − XiT γ∞ − Di XiT βc )
Ti Xi (Yiobs − XiT γ∞ − Di XiT βc )
n
i=1

−1  −1 Pn
P
00
+ n−1 ni=1 (1 − Ti )Xi e00i (0)
n
B
i=1 Ti Xi ei (1)
P
n−1 ni=1 Ti Xi e00i (1)
D

−1  −1 Pn
Pn
00 (0)
00
00 (0)} + n−1
X
e
−
e
n
B
i
i
i
i=1
i=1 Ti Xi {ei (1)
P
,
n−1 ni=1 Ti Xi e00i (1)
D
B
D

(A.18)
with (A, B, C, D) defined in (A.15) and {e00i (1), e00i (0)} defined in Theorem 7, and finally recognize
that the difference between the above two formulas has high order. Again we need only to find
the asymptotic distribution of ∆TSLS . The covariance of the second term on the right hand side of
(A.18) is (dropping the constant sum of Xi e00 (0))
 −1 Pn

00 (1) − e00 (0)}
n
T
X
{e
i
i
i
i
i=1
P
cov
n−1 ni=1 Ti Xi e00i (1)

1 n1 n0
S(Xε)
=
1
00 (1)} − S{Xe00 (0)} + S(Xε)]
2
[S{Xe
n n
2
11

1
00
2 [S{Xe (1)}


− S{Xe00 (0)} + S(Xε)]
,
S{Xe00 (1)}

where the off-diagonal term comes from the finite population covariance between X{e00 (1) − e00 (0)}
and Xe00 (1). Therefore, according to (A.18), the asymptotic covariance of ∆TSLS is the (2, 2) block
of the following matrix

−1
1 n1 n0 A B
·
C D
n2 n

S(Xε)
1
00 (1)} − S{Xe00 (0)} + S(Xε)]
[S{Xe
2

 
−T
− S{Xe00 (0)} + S(Xε)]
A B
·
,
S{Xe00 (1)}
C D

1
00
2 [S{Xe (1)}

which is
1 n1 n0 n −1 −1 −1
−1 −1 −1
−1 −1 −1 −1 T
−1 T
00
(p0 πc Sxx,c )S(Xε)(p0−1 πc−1 Sxx,c
) + (p−1
1 p0 πc Sxx,c )S{Xe (1)}(p1 p0 πc Sxx,c )
n2 n
o
−1 −1 −1 −1 T
−1 −1
00
00
−(p−1
0 πc Sxx,c )[S{Xe (1)} − S{Xe (0)} + S(Xε)](p1 p0 πc Sxx,c )


S{Xe00 (1)} S{Xe00 (0)} S(Xε)
= (πc Sxx,c )−1
+
−
(πc Sxx,c )−1 .
n1
n0
n
The asymptotic normality follows from the representation in (A.18) and the FPCLT.
Theorem 8: Decomposition of variation in non-compliance. The following proof uses two
facts: τa = τn = 0, and τ = πc τc .
Proof of Theorem 8. Write the total treatment effect variation as
Sτ τ

=

=

n

n

i=1
n
X

i=1

1X
1X 2
(τi − τ )2 =
τi − τ 2
n
n
1
n

I(Ui =c) τi2

i=1

−

πc2 τc2

= πc

n
1 X
I(Ui =c) τi2 − τc2
nc

!
+ πc (1 − πc )τc2 ,

i=1

the treatment effect variation explained by compliance status as
Sτ τ,U

=
=

X

πu (τu − τ )2 = πc (τc − πc τc )2 + πa (0 − πc τc )2 + πn (0 − πc τc )2

u=c,a,n

πc τc2 (1

− πc )2 + πc (πa + πn ) = πc (1 − πc )τc2 ,

and the subtotal treatment effect variation for compliers as
Sτ τ,c

n
n
1 X
1 X
2
I(Ui =c) (τi − τc ) =
I(Ui =c) τi2 − τc2 .
=
nc
nc
i=1

i=1

Therefore, the above three terms has the relationship Sτ τ = πc Sτ τ,c + Sτ τ,U .
The decomposition Sτ τ,c = Sδδ,c + Sεε,c follows immediately from the definition of βc .
12

Appendix B

More detailed comments

Appendices B.1–B.5 give more details of some technical issues and extensions mentioned in the
main text, and Appendix B.6 contains the proofs of the results in Appendix B.

Appendix B.1

Covariate adjustment to improve efficiency

In the main text, the role of covariates has been to model the treatment effect alone. In general,
we also want to use covariates to reduce sampling variability of βbRI , just as we can use covariates
to get more precise estimates of the average treatment effect. In particular, the goal is to more
bxt ∈ RK ; because these are the only random components in βbRI , if we estimate
precisely estimate S
them more precisely, we estimate βbRI more precisely as well. Let Wi ∈ RJ denote a vector of
pretreatment covariates without the intercept term. Because Xi and Wi have different roles in
estimation, they may also contain different sets of covariates, though, in practice, X is likely to be
a subset of W .
Following the covariate adjustment approach in survey sampling, we can obtain a model-assisted
estimator for β that uses W to reduce sampling variability. To see this, we need several definitions.
P
P
bww,t
Define W = n−1 ni=1 Wi and Sww = n−1 ni=1 Wi WiT , with det(Sww ) > 0; define W t and S
b t ∈ RJ×K as the regression
as the sample mean and covariance of W under treatment arm t; define B
coefficient of Y obs X on W for treatment arm t:
(
)
n
X
1
bt = S
b−1
B
I(Ti =t) Wi (Yiobs Xi )T .
ww,t
nt
i=1

The model-assisted estimator for Sxt is then
bw = S
bxt − B
b T (W̄t − W̄ ),
S
xt
t

(t = 0, 1).

As a result, we can improve the randomization-based estimator by
w
−1 bw
bw ).
βbRI
= Sxx
(Sx1 − S
x0
w is consistent for β with asymptotic covariance
Theorem A.1. The model-assisted estimator βbRI


S{E(0)} S(∆)
−1 S{E(1)}
−1
Sxx
+
−
Sxx
,
n1
n0
n

where Ei (t) = Yi (t)Xi − BtT (Wi − W̄ ) is the residual term and ∆i = Ei (1) − Ei (0).
13

w uses covariates both to estimate treatment effect variation and to reduce
The estimator, βbRI

sampling variability. Asymptotically, as long as W is predictive of the marginal potential outcomes,
the model-assisted estimator will improve precision over the unassisted estimators.

Appendix B.2

Fisherian exact inference

When εi = 0 for all i, we can obtain exact inference for β based on the Fisher randomization test
(Rubin, 1980; Rosenbaum, 2002; Ding et al., 2016). With a known β, the null hypothesis
H0 (β) : Yi (1) − Yi (0) = XiT β for all i

(A.19)

is sharp in the sense of allowing for full imputation of all missing potential outcomes based on
the observed data. We can perform randomization test using any sensible test statistic measuring
the deviation from the null hypothesis H0 (β), for example, the test statistic t(T , Y obs ; β) can be
the difference-in-means, difference-in-medians or the Kolmogorov–Smirnov statistics comparing two
samples {Yiobs − XiT β : Ti = 1, i = 1, . . . , n} and {Yiobs : Ti = 0, i = 1, . . . , n}. Then we can obtain
a (1 − α) level confidence region for β by inverting a sequence of randomization tests:
CRα = {β : Randomization test fails to reject H0 (β) at significance level α}.
The confidence region CRα is exact regardless of the sample size, and it is valid for general designs
of experiments if we use the corresponding assignment mechanism to simulate the null distribution
of the test statistic. Due to the duality between testing and interval estimation, we reject H0 (X)
with β1 = 0 in Section 3.3 if CRα ∩ {β : β1 = 0} is an empty set, which controls the type one error
rate by α.

Appendix B.3

A Variance Ratio Test

Raudenbush and Bloom (2015) have noticed that if the variance of the treatment potential outcome
is smaller than the control potential outcome, then the correlation between the individual treatment
effect and the control potential outcome is negative. This statement does not involve any covariates,
but it can be generalized to incorporate systematic and idiosyncratic treatment effect variation.
Below we give a finite population version of their result.
Theorem A.2. If the finite population variance of {Yi (1) − Xi0 β}ni=1 is smaller than {Yi (0)}ni=1 ,
then the idiosyncratic treatment effect variation, {εi }ni=1 , is negatively correlated with the control
potential outcomes.
14

Because the condition in Theorem A.2 depends only on the marginal distributions of the potential outcomes, we propose a formal variance ratio test of it using the observed data, which is a
generalization of a similar theorem in Ding et al. (2016):
Theorem A.3. The variance ratio test with rejection region
log s21 − log s20
p
< Φ−1 (α),
(b
κ1 − 1)/n1 + (b
κ0 − 1)/n0
has size at least as large as α, where s21 and κ
b1 are the sample variance and kurtosis of {Yiobs −
XiT βbRI : Ti = 1, i = 1, . . . , n}, and s20 and κ
b0 are the sample variance and kurtosis of {Yiobs : Ti =
0, i = 1, . . . , n}, and Φ−1 (α) is the α-th quantile of the standard normal distribution.
For finite population inference, the above test in Theorem A.3 is generally conservative, but for
superpopulation inference, it is asymptotically exact.
Note that Raudenbush and Bloom (2015) and Theorem A.2 are only about detecting a negative
association. Unfortunately, there is no testable condition for a positive association.

Appendix B.4

More on noncompliance: estimating the bounds of the R2 s

The component Sτ τ,U and and the probability πc are directly identifiable according to previous
discussion. Furthermore, Sδδ,c is also identifiable according to the following result.
Corollary A.1. Sδδ,c can be expressed as the expectation of the following quantity:
( n
)
n
n
X
X
1
1
1 1X
(δi − τc )2 −
Ti (1 − Di )(δi − τc )2 −
(1 − Ti )Di (δi − τc )2 .
πc n
n1
n0
i=1

i=1

i=1

T

Because πc , δi = Xi βc and τc can be estimated by a plug-in approach, Sδδ,c can also be estimated
from the observed data.
In the ITT case, estimation of the residual distributions are straightforward. In the noncompliance case, however, we need more discussion about the estimation of F1c (y) and F0c (y), because Ui
is a latent variable. To avoid notational clatter, we assume that γc1 and γc0 are known; in practice
bc1,RI and γ
bc0,RI , and the consistency
we can replace them by the randomization-based estimators γ
of the final estimator will not be affected. Recall the potential residuals e0i (1) and e0i (0) defined in
(17), and its observed value e0i = Ti e0i (1) + (1 − Ti )e0i (0). We define the following quantities
P
P
Fb10 (y) = n11 ni=1 Ti (1 − Di )I(e0i ≤y) ,
Fb11 (y) = n11 ni=1 Ti Di I(e0i ≤y) ,
P
P
(A.20)
Fb01 (y) = n10 ni=1 (1 − Ti )Di I(e0i ≤y) ,
Fb00 (y) = n10 ni=1 (1 − Ti )(1 − Di )I(e0i ≤y) .
Similar to Corollary 3, we have the following results.
15

Corollary A.2. For any y,
E{Fb11 (y) − Fb01 (y)} = πc F1c (y),

E{Fb00 (y) − Fb10 (y)} = πc F0c (y).

Therefore, we can estimate F1c (y) by {Fb11 (y) − Fb01 (y)}/b
πc , and estimate F0c (y) by {Fb00 (y) −
Fb10 (y)}/b
πc . As we mentioned before, in practice, we use eb0i instead of e0i in the formulas in (A.20).

Appendix B.5

Proofs of the theorems and corollaries in Appendix B

Proof of Theorem A.1. The population-level OLS regression matrix of Y (t)X onto W is
( n
)
1X
−1
T
Bt = Sww
Wi {Yi (t)Xi } ∈ RJ×K .
n
i=1

ew − S
ew ). According to the same argument
bxt + B T (W̄ − W̄t ) and βew = S −1 (S
ew = S
Define S
t
xt
xx
x1
x0
RI
w have the same asymptotic covariance, and in the following we need only to
as (A.4), βbRI and βeRI
w . Because
discuss the covariance of βeRI

ew − S
ew =
S
x1
x0

n
n

1 X 
1 X
Ti Yi (1)Xi + B1T (W̄ − Wi ) −
(1 − Ti ) Yi (0)Xi + B0T (W̄ − Wi )
n1
n0
i=1

=

i=1

n
n
1 X
1 X
Ti Ei (1) −
(1 − Ti )Ei (0)
n1
n0
i=1

i=1

can be represented as the difference between the sample means of Ei (1) and Ei (0), applying Theorem 2 we can obtain its covariance:


ew − S
ew = S{E(1)} + S{E(0)} − S{∆} ,
cov S
x1
x0
n1
n0
n
which completes the proof.
Proof of Theorem A.2. For simplicity, we abuse the variance and covariance notation for finite
Pn
2
T
population. For example, var{Y (0)} =
i=1 {Yi (0) − Ȳ (0)} /(n − 1). If var{Y (1) − X β} ≤
var{Y (0)}, then var{Y (0) + ε} ≤ var{Y (0)}. Expanding the left hand side,
var{Y (0)} + var{ε} + 2cov{Y (0), ε} ≤ var{Y (0)},
which implies 2cov{Y (0), ε} ≤ −var{ε} < 0.

16

Although it is straightforward to prove the conclusion for super population inference of Theorem
A.3 by using Ding et al. (2016, Theorem 2, Supplementary Material) and Slutsky’s Theorem, it
is less obvious to prove the conclusion for finite population inference. To simplify the proof, we
first prove the following lemma. Let (c1 , · · · , cn )T and (d1 , . . . , dn )T be two vectors of nonnegative
constants with the same mean m > 0 but different variances Scc and Sdd . The difference vector
(c1 − d1 , . . . , cn − dn )T has mean zero and variance Sc−d,c−d . Let
n
1 X
θbc =
Ti ci ,
n1
i=1

n
1 X
θbd =
(1 − Ti )di
n0
i=1

be two sample means of the treatment and control group, respectively.
Lemma A.5. Under the regularity conditions for the FPCLT, log θbc − log θbd has asymptotic mean
zero and variance
1
m2



Scc Sdd Sc−d,c−d
+
−
n1
n0
n


.

(A.21)

Proof of Lemma A.5. According to the FPCLT, we have the following joint asymptotic normality
of θbc and θbd :
θbc
θbd

!



  
Pn
Vcc Vcd
Ti ci
m
n−1
a
1
i=1
Pn
,
,
∼N
Vcd Vdd
m
n−1
0
i=1 (1 − Ti )di


=

where
n0
Scc ,
n1 n

n1
1
Sdd , Vcd = − (Scc + Sdd − Sc−d,c−d ).
n0 n
2n
b
b
Applying Taylor expansion at m, we have log θc −log θd = {(θbc −m)−(θbd −m)}/m+oP (n−1/2 ), which,
Vcc =

Vdd =

coupled with Neyman (1923)’s variance formula, gives the asymptotic variance of log θbc − log θbd in
(A.21).
Proof of Theorem A.3. First, as a direct consequence of Lemma A.5, the finite sample variance is
always larger than the super population variance, unless Sc−d,c−d = 0. Therefore, we need only to
show that the test in Theorem A.3 is asymptotically exact for super population inference, and the
asymptotic size of the test is no larger than α for finite population inference.
Second, replacing β by its consistent estimator βbRI does not affect the asymptotic distribution of
the test statistic, due to Slutsky’s Theorem. For simplicity, we treat β as known in our asymptotic
analysis.
With the two ingredients above, Theorem A.3 follows directly from the variance ratio test in
Ding et al. (2016, Theorem 2, Supplementary Material).
17

Proof of Corollary A.1. The conclusion follows from
(
)
(
)
n
n
n
1 X
1 X
1X
2
2
E
Ti (1 − Di )(δi − τc )
= E
Ti I(Ui =n) (δi − τc )
=
I(Ui =n) (δi − τc )2 ,
n1
n1
n
i=1
i=1
i=1
(
)
(
)
n
n
n
X
X
1
1
1X
2
2
E
(1 − Ti )Di (δi − τc )
= E
(1 − Ti )I(Ui =a) (δi − τc )
=
I(Ui =a) (δi − τc )2 .
n0
n0
n
i=1

i=1

i=1

Proof of Corollary A.2. We rewrite
Fb11 (y) =

n
n
1 X
1 X
Ti I(Ui =c) I{ei (1)≤y} +
Ti I(Ui =a) I{ei (1)≤y} ,
n1
n1

Fb10 (y) =

n
1 X
Ti I(Ui =n) I{ei (1)≤y} ,
n1

Fb01 (y) =

n
1 X
(1 − Ti )I(Ui =a) I{ei (0)≤y} ,
n0

Fb00 (y) =

n
n
1 X
1 X
(1 − Ti )I(Ui =c) I{ei (0)≤y} +
(1 − Ti )I(Ui =n) I{ei (0)≤y} .
n0
n0

i=1

i=1

i=1

i=1

i=1

i=1

In the above formulas, the random components are the Ti ’s, and therefore, the corollary follows
from Lemma A.1 and the linearity of expectations.

18

