Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.
Unifying Design-based Inference: On Bounding and Estimating the Variance of any Linear Estimator in any Experimental Design WORKING PAPER 1 OF 4
Introduction
This paper provides a design-based framework for variance (bound) estimation in experimental analysis. Results are applicable to virtually any combination of experimental design, linear estimator (e.g., difference-in-means, OLS, WLS) and variance bound, allowing for unified treatment and a basis for systematic study and compairison of designs using matrix spectral analysis. A proposed variance estimator reproduces Eicker-Huber-White (aka. “robust", “heteroskedastic consistent", “sandwich", “White", “Huber-White", “HC", etc.) standard errors and “cluster-robust" standard errors as special cases. While past work has shown algebraic equivalences between design-based and the so-called “robust" standard errors under some designs, this paper motivates them for a wide array of design-estimator-bound triplets. In so doing, it provides a clearer and more general motivation for “robust" variance estimators.
The Neyman Causal Model
Consider a randomized experiment with $k$ treatment arms. The Neyman causal model (NCM) assumes that the units in the experimental study represent a finite population of size $n$. For a given outcome measure, call it $y$, each unit, $i$, responds with one of $k$ possible values in $\{y_{1i}, y_{2i}, ..., y_{ki} \}$, depending on their treatment assignment. The possible responses are referred to as the potential outcomes. In the NCM these values are considered (nonrandom) constants, which stands in contrast to other, more common, formulations where potential outcomes are assumed to be sampled from some (possibly nonparametric) distribution.
The only random element in the NCM is the treatment assignment indicators $\{R_{1i}, R_{2i},...,R_{ki} \}$, and they determine which potential outcome will be observed by the researcher. Since a unit can only be assigned to one arm of the experiment, only one of the indicators will realize a value of one, and the rest will be zero, such that $R_{1i}+R_{2i}+...+R_{ki}=1$ for all $i$.
A standard representation of the observed outcome for the $i^{th}$ unit under the NCM would be,
align*[align* omitted — 68 chars of source]
which is itself random, due to the assignment indicators. For each unit, the observed data can then be represented as $\{Y^{obs}_i, R_{1i}, R_{2i},..., R_{ki}, x_i\}_{\forall i}$, where $x_i$ is an additional vector of $l$ covariates. Like the potential outcomes, $x_i$ is nonrandom, but unlike the potential outcomes the same value is observed irrespective of the assignment.
Ideally, we would like to know, for a given individual, $i$, the difference between responses under various arms, called a treatment effect. It is clear from the definition of $Y_i^{obs}$, however, that individual treatment effects are not observable since only one of the potential outcomes can be observed for an individual, a problem known as fundamental problem of causal inference holland. As a result, researchers often try to estimate averages of across the units in study.
Example (Treatment/Control Experiment) : In an experiment with a control group (arm 0) and a treatment group (arm 1) the individual-level treatment effect, $y_{1i}-y_{0i}$, but this is not identified, so a researcher might try to estimate the average treatment effect $n^{-1} \sum_i \left(y_{1i}-y_{0i}\right)$.{ $\triangle$}
Example ($2 \times 2$ Factorial Experiment) : Consider a 2$\times$2 factorial design with treatments A and B. Units in arm 1 are controls (no treatments), units in arm 2 are given treatment A only, units in arm 3 are given B only, and units in arm 4 are given both A and B. Similar to the treatment/control example, one could contrast the mean of an arm with a single treatment against the control mean, e.g., the average effect of A compared to no treatments, $n^{-1} \sum_i \left(y_{2i}-y_{1i}\right)$. Another quantity of interest might be an average marginal causal effect (AMCE), e.g., the effect of A marginalizing over the levels of B, $n^{-1} \sum_i \frac{1}{2} \left(y_{2i}-y_{1i}+y_{4i}-y_{3i}\right)$. Another example might be an omnibus test based on the contrast $n^{-1} \sum_i \left(y_{2i}/3+y_{3i}/3+y_{4i}/3-y_{1i}\right)$.{ $\triangle$}
Target quantities such as local average treatment effects or conditional average treatment effects might also be considered in this framework, but the primary focus of this paper is variance estimation for linear estimators for virtually any design.
Suffice to say that developing variance estimators before considering point estimation is appealing, if somewhat counter-intuitive, for two reasons. On the one hand, asymptotic analysis for point estimators can be made easier by having first established general variance expressions (for all linear estimators and virtually any design). On the other hand,
a general framework for variance (bound) estimation can be developed even while a particular estimation target has yet to be defined, and even if an “estimator" does not estimate anything of interest, it's variance can still be studied.
Notation
To simplify notation, let $y_1$, $y_2$,...,$y_k$ represent length $n$ vectors of potential outcomes associated with each of the arms, with the $i^{th}$ element of each corresponding to the $i^{th}$ unit. Next, stack these vectors to create
align*[align* omitted — 95 chars of source]
which is a column vector and has length $kn$ containing all $k$ potential outcomes for all $n$ units.
Next, if we let $1_{\scriptscriptstyle n}$ be a $n$-length vector of ones, then a $kn \times k$ intercept matrix can be defined as,
align*[align* omitted — 194 chars of source]
which, for example, allows us to express a k-length vector of the means of each arm as $\frac{1}{n} \mathds{1}' y$, or, equivalently, $\left(\mathds{1}' \mathds{1}\right)^{-1} \mathds{1}' y$. Next, define $c$ as the contrast vector, of length $k$, such that $c'\left(\mathds{1}' \mathds{1}\right)^{-1} \mathds{1}' y$ gives contrasts between potential outcome means for the various arms.
Example (Treatment/Control Experiment, continued) :
With two arms, control (arm 1) and treatment (arm 2), define $c=\left(-1 \hspace{2mm} 1\right)'$. Then the average treatment effect is simply $c' \left(\mathds{1}' \mathds{1}\right)^{-1} \mathds{1}' y= n^{-1} \sum_i \left(y_{2i}-y_{1i}\right)$.{ $\triangle$}
Example ($2 \times 2$ Factorial Experiment, continued) : In a four-arm experiment, if $c=(-1\hspace{2mm} 1 \hspace{2mm} 0 \hspace{2mm} 0)'$ then $ c' \left(\mathds{1}' \mathds{1}\right)^{-1} \mathds{1}' y= n^{-1} \sum_i \left(y_{2i}-y_{1i}\right) $ is the avearge difference between the first two arms. Alternatively, if the researcher chooses $c=(-\frac{1}{2}\hspace{2mm} \frac{1}{2} \hspace{2mm} -\frac{1}{2} \hspace{2mm} \frac{1}{2} )'$ then $c' \left(\mathds{1}' \mathds{1}\right)^{-1} \mathds{1}' y = n^{-1} \sum_i \frac{1}{2} \left(y_{2i}-y_{1i} + y_{4i}-y_{3i} \right)$.
{ $\triangle$}
Next define an $n \times n$ diagonal matrix that has all $n$ assignment indicators for treatment arm 1 on the diagonal,
align*[align* omitted — 174 chars of source]
and define $\mathbf{R}_2$, $\mathbf{R}_3$, $\hdots$, $\mathbf{R}_k$ analogously. Arrange these matrices to create the diagonal $kn \times kn$ matrix
align*[align* omitted — 156 chars of source]
and note the a $kn\times kn$ diagonal matrix of assignment probabilities can be written as $\boldsymbol{\pi}:=\text{\textnormal{E}}[\mathbf{R}]$, with the first $n$ diagonal elements representing probabilities of assignment to arm 1, then the next $n$ diagonal elements are probabilities of assignment to arm 2 and so on.
In this alternative notation the researcher can be said to observe the assignment, $\mathbf{R}$, the observed vector of outcomes, $\mathbf{R} y$, and also a matrix of $l$ pre-treatment covariates, $\mathbf{x}$, which has size $n \times l$. In a randomized experiment $\boldsymbol{\pi}$ is also observed (known) in many cases. When intractable analytically, however, it might be estimated to arbitrary precision by repeating the original randomization until a target level of precision is achieved.
For covariate adjusted estimators, it is convenient to define the $kn \times (k+l) $ matrix,
align*[align* omitted — 240 chars of source]
which augments the intercept vector, $\mathds{1}$, with covariates.
remarkFor some cases, such adjusting for covariates separately by arm, it might be useful to define $\mathbb{x}$ with $\mathbf{x}$ matrices arranged along a block-diagonal. In that case, it is useful to stipulate that $\mathbf{x}$ have columns that sum to zero to avoid problems of coefficient interpretation lin, middleton18. This will be discussed further in paper 3 of 4.
Linear estimators
This paper focuses on the variance, bounding and variance bound estimation of the class of estimators that are linear in the observed outcome, $y$. This class includes everything from the difference-of-means, to the Horvitz-Thomposon estimator, to regression.
Note, however, that beyond presenting a general approach to variance bound estimation for the class of linear estimators, point estimation itself will be the focus of the third and fourth papers in the series. Questions such as consistency will be and causal identification will be considered then. For now, suffice it to be said that an estimator need not be consistent for any quantity of interest at all (causal or otherwise) in order to derive variance expressions for it.
Definition
definition[Linear Estimators]
Linear estimators are defined as having the form,
\begin{align}
\widehat{\delta}_c := & c' \mathbf{W} \mathbf{R} y,
\end{align}
where $\mathbf{W}$ a matrix with $kn$ columns and $k$ rows if it is an unadjusted estimator and $k+l$ rows if it is a covariate adjusted estimator. The length of the contrast vector, $c$, is equal to the number of rows in $\mathbf{W}$. The first $k$ entries of $c$ are the contrast values, followed by $l$ zeros in the case of covariate adjusted estimators.
Also, for convenience, define $\mathbf{w}$ to be $\mathbf{W}$ evaluated at $\mathbf{R}=\boldsymbol{\pi}$, i.e.,
align[align omitted — 85 chars of source]
definition[Horvitz-Thompson estimator]
The Horvitz-Thompson estimator written as in Definition (ref) with,
\begin{align*}
\mathbf{W} &= \mathbf{W}^{\scriptscriptstyle{HT}} :=\left(\mathds{1}'\mathds{1}\right)^{-1} \mathds{1}'\boldsymbol{\pi}^{-1} ,
\end{align*}
noting that $\mathbf{W}^{\scriptscriptstyle{\textnormal{HT}}}=\mathbf{w}^{{\scriptscriptstyle{\textnormal{HT}}}}$ since $\mathbf{W}^{\scriptscriptstyle{\textnormal{HT}}}$ is nonrandom.
definition[Contrast-of-means]
Contrast-of-means (e.g., difference-of-means) can be written as in Definition 2.1 with,
\begin{align*}
\mathbf{W} &= \mathbf{W}^{\scriptscriptstyle{CM}} := \left(\mathds{1}' \mathbf{R} \mathds{1} \right)^{-1}\mathds{1}' .
\end{align*}
definition[Hajek estimator]
The Hajek estimator can be written as Definition ((ref)) with,
\begin{align*}
\mathbf{W} &= \mathbf{W}^{\scriptscriptstyle{HJ}} := \left(\mathds{1}' \boldsymbol{\pi}^{-1} \mathbf{R} \mathds{1} \right)^{-1}\mathds{1}' \boldsymbol{\pi}^{-1} .
\end{align*}
definition[OLS estimator]
The OLS estimator can be written as Definition ((ref)) with,
\begin{align*}
\mathbf{W} &= \mathbf{W}^{\scriptscriptstyle{OLS}} := \left(\mathbb{x}' \mathbf{R} \mathbb{x} \right)^{-1}\mathbb{x} ' .
\end{align*}
definition[WLS estimators]
WLS estimators can be written as in Definition 2.1 with,
\begin{align*}
\mathbf{W} &= \mathbf{W}^{\scriptscriptstyle{WLS}} := \left(\mathbb{x}' \mathbf{m} \mathbf{R} \mathbb{x} \right)^{-1}\mathbb{x}' \mathbf{m} .
\end{align*}
remarkWLS is a class that includes OLS, Hajek and contrast-of-means (e.g., difference-of-means) as special cases. It is equivalent to OLS when $\mathbf{m}=\mathbf{i}_{kn}$ ($\mathbf{i}_{kn}$ is the identity matrix). If $\mathbf{m}=\mathbf{i}_{kn}$ and, in addition, $\mathbb{x}=\mathds{1}$, WLS is OLS without covariates, which is equivalent to the contrast-of-means (e.g., in the two-arm case, we call this the difference-of-means), underscoring Theorem 1 in freedman08a. If $\mathbf{m}=\boldsymbol{\pi}^{-1}$ and $\mathbb{x}=\mathds{1}$, then it is the Hajek estimator. The covariate adjusted WLS with $\mathbf{m}=\boldsymbol{\pi}^{-1}$ will be discussed further in paper 3 of 4, because it is algebraically equivalent to the generalized regression estimator introduced there.
First-order Taylor approximation
In this section, a general approach to obtaining asymptotically valid variance expressions for linear estimators is given using a first-order approximation of a Tyalor series. The method is often used when an exact, closed-form variance expression is not tractable, as may be the case with any number of linear estimators. Examination of the $W$ vectors defined above shows that, with the exception of Horvitz-Thompson, the estimators all had random denominators (i.e., inverted random matrices), making closed form variance expressions difficult.
The original estimator and its Taylor approximation are asymptotically equivalent (cite Pashley). As such, the original estimator “borrows" the closed-form variance expression given for the Taylor approximation, again justified given the asymptotic equivalence.
lemma[First-order Taylor approximation for linear estimators]
First, assume a linear estimator as defined in Definition 2.1. Then, let $\Big \{ \left. . \hspace{1mm} \right| _{\mathbf{R} = \boldsymbol{\pi}} \Big \}$ represent a function that evaluates the argument to the left of the vertical line at ${\mathbf{R} = \boldsymbol{\pi}}$. Similarly, let $\Big \{ \left. . \hspace{1mm} \right| _{\mathbf{R} = \boldsymbol{\pi}} \left(\mathbf{R}-\boldsymbol{\pi}\right)\Big \}$ evaluate its argument at ${\mathbf{R} = \boldsymbol{\pi}}$ and then multiply by $\left(\mathbf{R}-\boldsymbol{\pi}\right)$. Then from Taylor's theorem and the product rule, we have the first-order Taylor approximation, $\widehat{\delta} \approx \widehat{\delta}^{\scriptscriptstyle{\textnormal{T}}}$, with
\begin{align}
\widehat{\delta}^{\scriptscriptstyle{T}}_c
:= & \Big\{ \left. c' \mathbf{W} \mathbf{R} y \right|_{ \mathbf{R} =\boldsymbol{\pi}}\Big \}
+\Big\{ \left. c' \mathbf{W} \right|_{ \mathbf{R} =\boldsymbol{\pi}} \Big \}\bigg \{ \left. \frac{d}{d \mathbf{R}} \mathbf{R} \right|_{ \mathbf{R} =\boldsymbol{\pi}} \left(\mathbf{R}-\boldsymbol{\pi}\right) \bigg \} y \nonumber
\\ &+ \bigg \{ \left. \frac{d}{d \mathbf{R}} c' \mathbf{W} \right|_{ \mathbf{R} =\boldsymbol{\pi}} \left(\mathbf{R}-\boldsymbol{\pi}\right) \bigg\} \Big\{ \left. \mathbf{R} \right|_{ \mathbf{R} = \boldsymbol{\pi} } \Big \}y \nonumber
\\ =& a_c + c' \mathbf{w} \mathbf{R} y + \bigg \{ \left. \frac{d}{\text{d} \mathbf{R}} c' \mathbf{W} \right|_{ \mathbf{R} =\boldsymbol{\pi}} \mathbf{R} \bigg\} \boldsymbol{\pi} y
\end{align}
where
\begin{align*}
a_c =& - \bigg \{ \left. \frac{\text{d}}{\text{d} \mathbf{R}} c' \mathbf{W} \right|_{ \mathbf{R} =\boldsymbol{\pi}} \boldsymbol{\pi} \bigg\} \boldsymbol{\pi} y
\end{align*}
is a constant.
remarkAn expression for $a_c$ is given but it is not important for the purposes of variance approximations because the term is a constant. Recall that the purpose of deriving a first-order Taylor approximation, $\widehat{\delta}^{\scriptscriptstyle{\textnormal{T}}}_c$, is to identify a closed-form variance expression that might then be “borrowed" by the original linear estimator given in Definition (ref).
theoremFor a constant, $a_c$, and vector of constants, $z_c$, first-order Taylor approximations for linear estimators may be written as,
\begin{align*}
\widehat{\delta}^{\scriptscriptstyle{T}}_c
= & a_c
+ n 1_{\scriptscriptstyle {k}}' \mathbf{w}^{{\scriptscriptstyle{HT}}} \mathbf{R} z_c,
\end{align*}
where $z_c$ has the form $z_c = \boldsymbol{\pi} \textnormal{\footnotesize diag}\left( \mathbf{l} y\right) \mathbf{t}' c$ and where $(k \times kn)$ matrix $\mathbf{t}$ and $(kn \times kn)$ matrix $\mathbf{l}$ depend on the estimator. Hence, a first-order approximation of a Tyalor series using Taylor's theorem variance approximations will be expressed as the variance of a Horvitz-Thompson estimator of the ATE of $z_c$ with contrast vector $n 1_{\scriptscriptstyle {k}}$.
proofWith Equation ((ref)), it is easy to see that the Taylor linearized approximation has the form
\begin{align*}
\widehat{\delta}^{\scriptscriptstyle{T}}_c= & a_c +c' \mathbf{t} \mathbf{R} \mathbf{l} y
\end{align*}
where matrices $\mathbf{t}$ and $\mathbf{l}$ are $(k \times kn)$ and $(kn \times kn)$, respectively, and will depend on the estimator. Noting that $c' \mathbf{t}$ is a $(1 \times kn)$ vector, write
\begin{align*}
\widehat{\delta}^{\scriptscriptstyle{T}}_c = & a_c + 1_{\scriptscriptstyle {kn}}' diag\left( c' \mathbf{t} \right) \mathbf{R} \mathbf{l} y
\\ = & a_c + 1_{\scriptscriptstyle {kn}}' \mathbf{R} diag\left( c' \mathbf{t} \right) \mathbf{l} y
\\ = & a_c + 1_{\scriptscriptstyle {kn}}' \boldsymbol{\pi}^{-1} \mathbf{R} \boldsymbol{\pi} diag\left( c' \mathbf{t} \right) \mathbf{l} y
\\ = & a_c + n 1_{\scriptscriptstyle {k}}' (\mathds{1}' \mathds{1})^{-1} \mathds{1}' \boldsymbol{\pi}^{-1} \mathbf{R} \boldsymbol{\pi} diag\left( c' \mathbf{t} \right) \mathbf{l} y
\\ = & a_c + n 1_{\scriptscriptstyle {k}}' \mathbf{w}^{\scriptscriptstyle{\textnormal{HT}}} \mathbf{R} z_c
\end{align*}
where $z_c := \boldsymbol{\pi} \textnormal{\footnotesize diag}\left( \mathbf{l} y\right) \mathbf{t}' c$.
remarkThe result shows that first order Taylor approximations are Horvitz-Thompson estimators. This highlights the importance of studying Horvitz-Thompson variance in order to develop asymptotic variance expressions for linear estimators in general.
remarkThe constant vector $z_c$ is not directly observed. The next section will show that the plug-in principle provides a basis for asymptotically valid variance expressions.
{12pt}
table[table omitted — 3,197 chars of source]
{0pt}
Example (Weighted least squares) : Weighted least squares is a class that includes OLS ($\mathbf{m}=\mathbf{i}_{kn}$), contrast-of-means (e.g., difference of means, with $\mathbf{m}=\mathbf{i}_{kn}$ and $\mathbb{x}=\mathds{1}$) and the Hajek estimator ($\mathbf{m}=\boldsymbol{\pi}^{-1}$ and $\mathbb{x}=\mathds{1}$). To derive its Taylor approximation, first let $\mathbf{w}^{\scriptscriptstyle{\textnormal{WLS}}}= {\left. {\mathbf{W}^{{\scriptscriptstyle{\textnormal{WLS}}}}} \right| }_{\mathbf{R}=\boldsymbol{\pi}}= \left( \mathbb{x}' \mathbf{m} \boldsymbol{\pi} \mathbb{x}\right)^{-1} \mathbb{x}' \mathbf{m} $, and note that by the rules of matrix differentiation the third term in Equation ((ref)) is
align*[align* omitted — 1,133 chars of source]
Therefore, Equation ((ref)) made specific to WLS is
align*[align* omitted — 1,308 chars of source]
where $z_c^{{\scriptscriptstyle{\textnormal{WLS}}}} = \boldsymbol{\pi} \textnormal{\footnotesize diag}\left(y- \mathbb{x} b^{{\scriptscriptstyle{\textnormal{WLS}}}}\right) {\mathbf{w}^{\scriptscriptstyle{\textnormal{WLS}}}}' c $ is recognizable in the form given in Theorem (ref) with $\textbf{t}^{\scriptscriptstyle{\textnormal{WLS}}}=\mathbf{w}^{{\scriptscriptstyle{\textnormal{WLS}}}}$ and $\textbf{l}^{{\scriptscriptstyle{\textnormal{WLS}}}} = \mathbf{i}_{kn} -\mathbb{x} \left( \mathbb{x}' \mathbf{m} \boldsymbol{\pi} \mathbb{x} \right)^{-1} \mathbb{x}' \mathbf{m} \boldsymbol{\pi} $ is a “residual maker" matrix.
{ $\triangle$}
comment\begin{table}[H]
\begin{center}
\caption{Examples of linear estimators. $W$ is as defined in Definition 2.1, $z$ is as defined in Theorem (ref).}
\begin{tabular}{l | c | c | c }
& abbr. & $W$ & $z'$ \\
\hline Horvitz-Thompson & ${\scriptscriptstyle{\textnormal{HT}}}$ & $n^{-1} 1_{\scriptscriptstyle {kn}}' \boldsymbol{\pi}^{-1} $ & $y'$
\\
Difference-of-means & ${\scriptscriptstyle{\textnormal{CM}}}$ & $\big\{ \left(\mathds{1}' \mathbf{R} \mathds{1}\right)^{-1}\mathds{1}' \big\}_{\scriptscriptstyle (2,.)}$ & $ n \big\{ \mathbf{w}^{\scriptscriptstyle{\textnormal{CM}}} \big\}_{\scriptscriptstyle (2,.)} \textnormal{\footnotesize diag}(y - \mathds{1} \mathbf{w}^{\scriptscriptstyle{\textnormal{CM}}} \boldsymbol{\pi} y) \boldsymbol{\pi} $
\\
Hajek & ${\scriptscriptstyle{\textnormal{HJ}}}$ & $\big\{\left(\mathds{1}' \boldsymbol{\pi}^{-1} \mathbf{R} \mathds{1}\right)^{-1}\mathds{1}' \boldsymbol{\pi}^{-1} \big\}_{\scriptscriptstyle (2,.)}$ & $y - \mathds{1}\left( \mathds{1}' \mathds{1}\right)^{-1} \mathds{1}'y$
\\
OLS & ${\scriptscriptstyle{\textnormal{OLS}}}$ & $\big\{ \left(\mathbb{x}' \mathbf{R} \mathbb{x} \right)^{-1}\mathbb{x}' \big\}_{\scriptscriptstyle (2,.)}$ & $\textnormal{\footnotesize diag}\left(y- \mathbb{x} b^{{\scriptscriptstyle{\textnormal{OLS}}}}\right) \big\{ n \left( \mathbb{x}' \boldsymbol{\pi} \mathbb{x}\right)^{-1} \mathbb{x}' \boldsymbol{\pi} \big\}_{\scriptscriptstyle (2,.)}'$
\\
WLS ($\mathbf{m}=\boldsymbol{\pi}^{-1}$) & ${\scriptscriptstyle{\textnormal{WLS}}}$ & $\big\{\left(\mathbb{x}' \boldsymbol{\pi}^{-1} \mathbf{R} \mathbb{x} \right)^{-1}\mathbb{x}' \boldsymbol{\pi}^{-1} \big\}_{\scriptscriptstyle (2,.)}$ & $\textnormal{\footnotesize diag}\left(y- \mathbb{x} b^{{\scriptscriptstyle{\textnormal{WLS}}}}\right) \big\{ n \left( \mathbb{x}' \mathbb{x}\right)^{-1} \mathbb{x}' \big\}_{\scriptscriptstyle (2,.)}'$
\end{tabular} \end{center}
\end{table}
Variance
Now that the importance of Horvitz-Thompson estimators for asymptotic variance expressions for the entire class of linear estimators (which includes, for example, OLS, WLS, Hajek, and difference-of-means) has been established, this section will give the variance of HT estimators and first-order approximates of linear estimators.
Throughout, we will make use of the $kn \times kn$ “first order design matrix", which will allow for easy comparison of designs using spectral analysis.
definitionThe “first-order design matrix" is a variance-covariance matrix of inverse-probability weighted treatment assignments, written,
\begin{align}
\mathbf{d}:=&V \left( 1_{\scriptscriptstyle {kn}}'\boldsymbol{\pi}^{-1} \mathbf{R} \right)
\\ = & \left(E\left[\mathbf{R} 1_{\scriptscriptstyle {kn}} 1_{\scriptscriptstyle {kn}}' \mathbf{R} \right] - \boldsymbol{\pi} 1_{\scriptscriptstyle {kn}} 1_{\scriptscriptstyle {kn}}' \boldsymbol{\pi} \right) / \left( \boldsymbol{\pi} 1_{\scriptscriptstyle {kn}} 1_{\scriptscriptstyle {kn}}' \boldsymbol{\pi} \right), \nonumber
\end{align}
where “/" represents elementwise division.
theorem[Horvitz-Thompson Variance] An exact expression for the variance of Horvitz-Thompson estimators is given by
\begin{align*}
V\left(\widehat{\delta}^{{\scriptscriptstyle{HT}}}_c \right) = { z^{{\scriptscriptstyle{HT}}}_c}'\mathbf{d} { z^{{\scriptscriptstyle{HT}}}_c},
\end{align*}
where ${z_c^{{\scriptscriptstyle{\textnormal{HT}}}}}' = c' \left(\mathds{1}' \mathds{1} \right)^{-1} \mathds{1}' \textnormal{\footnotesize diag}\left(y\right)$.
proofUsing the identity $y=\textnormal{\footnotesize diag}\left(y\right) 1_{\scriptscriptstyle {kn}}$, the Horvitz-Thompson estimator can be written
\begin{align*}
\widehat{\delta}^{{\scriptscriptstyle{HT}}}
= & c' \left(\mathds{1}' \mathds{1} \right)^{-1} \mathds{1}' \boldsymbol{\pi}^{-1} \mathbf{R} diag\left(y\right) 1_{\scriptscriptstyle {kn}}
\\ = & c' \left(\mathds{1}' \mathds{1} \right)^{-1} \mathds{1}' diag\left(y\right) \boldsymbol{\pi}^{-1} \mathbf{R} 1_{\scriptscriptstyle {kn}}.
\\ = & { z^{{\scriptscriptstyle{HT}}}_c}' \boldsymbol{\pi}^{-1} \mathbf{R} 1_{\scriptscriptstyle {kn}}.
\end{align*}
where ${z_c^{{\scriptscriptstyle{\textnormal{HT}}}}}' = c' \left(\mathds{1}' \mathds{1} \right)^{-1} \mathds{1}' \textnormal{\footnotesize diag}\left(y\right)$.
So the variance can be written,
\begin{align*}
V\left(\widehat{\delta}^{{\scriptscriptstyle{\textnormal{T}}}}_c \right) & = { z^{{\scriptscriptstyle{\textnormal{HT}}}}_c}'\text{\textnormal{V}} \left( 1_{\scriptscriptstyle {kn}}'\boldsymbol{\pi}^{-1} \mathbf{R} \right) { z^{{\scriptscriptstyle{\textnormal{HT}}}}_c}
\\ & = { z^{{\scriptscriptstyle{\textnormal{HT}}}}_c}'\mathbf{d} { z^{{\scriptscriptstyle{\textnormal{HT}}}}_c}
\end{align*}
theorem[Variance of first-order Taylor approximations] The variance of first-order Taylor approximations of linear estimators can be written as,
\begin{align}
V\left(\widehat{\delta}^{{\scriptscriptstyle{T}}}_c \right) = {z_c^{\scriptscriptstyle{T}}}' \mathbf{d} z_c^{\scriptscriptstyle{T}},
\end{align}
with examples of $z^{\scriptscriptstyle{\textnormal{T}}}_c$ given in Table (ref).
proofBy Theorem (ref), linear approximations are Horvitz-Thompson estimators of a vector $z_c$ and contrast vector $n 1_{\scriptscriptstyle {k}}$. Now, $ n 1_{\scriptscriptstyle {k}}' \left( \mathds{1}' \mathds{1} \right)^{-1} \mathds{1}' \textnormal{\footnotesize diag}\left(z^{\scriptscriptstyle{\textnormal{T}}}_c\right)= 1_{\scriptscriptstyle {kn}}' \textnormal{\footnotesize diag}\left(z^{\scriptscriptstyle{\textnormal{T}}}_c\right) ={z^{\scriptscriptstyle{\textnormal{T}}}_c} $. Therefore, using Theorem (ref),
\begin{align*}
V\left(\widehat{\delta}^{{\scriptscriptstyle{T}}}_c \right) & = n 1_{\scriptscriptstyle {k}}' \left(\mathds{1}' \mathds{1} \right)^{-1} \mathds{1}' diag\left(z^{\scriptscriptstyle{T}}_c\right) \mathbf{d} diag\left(z^{\scriptscriptstyle{\textnormal{T}}}_c\right) \mathds{1} \left(\mathds{1}' \mathds{1} \right)^{-1} 1_{\scriptscriptstyle {k}} n
\\ & = {z^{\scriptscriptstyle{\textnormal{T}}}_c}' \mathbf{d} {z^{\scriptscriptstyle{\textnormal{T}}}_c}
\end{align*}
remarkEquation ((ref)) is an exact variance expression, however, they are not identified. The next subsection introduces the necessary concept of variance bounding.
remarkThe design matrix, $\mathbf{d}$, will provide useful device in the study the best designs for outcomes with different characteristics as the next example will show.
Example (Comparing complete randomization and paired randomization) : Consider a treatment/control (two-arm) experiment that is pair-randomized. A pair-randomized design is a special case of a block-randomized (i.e., stratified) design where blocks have size 2. In each pair/block, one unit is assigned to treatment and the other in control with equal (.5) probability. Across blocks, assignments are independent.
When $n=4$ (and assuming w.l.o.g. that the data are sorted by pair), the design matrix is
align*[align* omitted — 510 chars of source]
and note that empty cells are 0. The design matrix for complete randomization (where 2 of 4 are randomly assigned to treatment) is
align*[align* omitted — 1,699 chars of source]
Eigendecomposition of $\mathbf{d}^{cr}-\mathbf{d}^{pr}$ gives eigenvalues $2.67, 0, 0, 0, 0, 0, -1.33,$ and $-1.33$ corresponding eigenvectors in Table (ref). The eigenvectors associated with nonzero eigenvalues provide insight into the subspace in $ {\mathds {R}}^{\scriptscriptstyle 2n}$ where one design may be preferable to another, for example, when the estimator is difference-in-means (which is equivalent to both Horvitz-Thopson and Hajek for these designs).
table[table omitted — 893 chars of source]
In this example, assuming the contrast matrix is $c=\left(-1,1\right)'$, and examining the first eigenvector with eigenvalue 2.67, one can conclude that if the outcomes for the four units given in Table (ref), then the difference-of-means would be much less precise under the completely randomized design. So, the eigenvector in a sense represents a “best-case" (normed) potential outcome vector for paired randomization. Inspection of the outcomes themselves confirms the intuition that pair randomization is better than complete randomization when units are homogenous within pairs.
table[table omitted — 384 chars of source]
Next, considering the two eigenvectors associated with the eigenvalue -1.33, we see the implied potential outcomes in Table (ref) give potential outcomes for which complete randomization is preferable. Note that either of the two sets is a “worst-case" scenario for paired randomization, as is any set of potential outcomes that can be generated by linear combinations of the two eigenvectors. Inspection of these outcomes is consistent with the observation that complete randomization can be better than paired randomization when paired units are maximally heterogeneous.
table[table omitted — 518 chars of source]
{ $\triangle$}
remarkThe example illustrates a relatively effortless method of identifying key insights about arbitrary designs through spectral analyses of first-order design matrices. In the example, the observation that pair randomization can hurt precision when units are not homogeneous within pairs is not new. However, this approach to comparing designs is perfectly general and can be applied to virtually any designs.
Variance bounds
In spite of an exact expression for first-order Taylor approximations in Equation ((ref)), the quantity is never identified because not all terms in the quadratic can be observed. Even if the elements of $\mathbf{R} z_c$ were observed directly (which is the case for $\mathbf{R} z_c^{\scriptscriptstyle{\textnormal{HT}}}$ but none of the other examples in Table 1), some pairs of potential outcomes can never be jointly observed. For example, for a given unit, only one of two (or more) potential outcomes can be observed, a problem is referred to as the “fundamental problem of causal inference" holland. Other design features, such as clustering or pair randomization, can also render various combinations of potential outcomes unobservable.
Starting with Neyman (1923) one proposed solution to unidentified variance has been to estimate a {\it variance bound}, i.e., a quantity that is provably greater than the variance, but which is identified. It should be understood that while the term variance estimation is often used as a shorthand in the literature, it is not, in general, an accurate phrase. Variance bound estimation is a more precise so it will be used here.
definition[Variance bound matrix]
Let $\tilde{\mathbf{d}}$ be an arbitrary $kn \times kn$ matrix and let be $z$ an arbitrary vector with length $kn$. Then $\tilde{ \mathbf{d} }$ is a variance bound matrix (or bounding matrix) for $\mathbf{d}$ if, for all $z\in\mathds{R}^{kn}$, $z'\mathbf{d} z \leq z'\tilde{\mathbf{d}}z$.
lemma$\tilde{ \mathbf{d} }$ is a bounding matrix $\mathbf{d} $ if and only if matrix $\tilde{ \mathbf{d} }-\mathbf{d}$ is positive semi-definite.
proofBy the definition of a bound, $ z'\tilde{\mathbf{d}}z - z'\mathbf{d} z \geq 0$ for all $z\in \mathds{R}^{kn}$. This implies that $z'(\tilde{\mathbf{d}}-\mathbf{d} ) z \geq 0$, i.e., that $ \tilde{\mathbf{d}}-\mathbf{d} $ is positive semi-definite.
definition[Identified variance bound]
Let $\tilde{\mathbf{d}}$ be bounding matrix for $\mathbf{d}$. It gives an identified variance bound if
\begin{align*}
I (\mathbf{d}=-1)
\circ I(\tilde{\mathbf{d}}= 0)=\textnormal{I} (\mathbf{d}=-1)
\end{align*}
where $\circ$ is element-wise multiplication, $\text{\textnormal{I}}\left(\mathbf{d}=-1\right)$ is an indicator function returning an $kn \times kn$ matrix of ones and zeros indicating whether each element of $\mathbf{d}$ is equal to $-1$ (an indication that the associated term in the variance quadratic is impossible to observe), and $\text{\textnormal{I}} (\tilde{ \mathbf{d} } =0 )$ is, similarly, an indicator function returning an $kn \times kn$ matrix of ones and zeros indicating the location of zeros in $\tilde{ \mathbf{d} }$.
Generalizing Neyman's variance bound
This section proposes a generalization of Neyman's (1923) variance bound. Let matrix $\mathbf{d}$ be partitioned into $k^2$ partitions of size $n \times n$. Then for $r,s \in \{1, 2, ..., k\}$, let the $\mathbf{d}_{rs}$ be the $(r,s)^{th}$ partition, having dimension $n \times n$. Also, let $c_r$ be the $r^{th}$ element of the length-$k$ contrast vector, $c$.
Then the following bounding method produces an identified bound for experiments when partitions $\textnormal{I}(\tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{N}}}_{rr}==-1)=0_{\scriptscriptstyle n\times n}$, i.e., there are no $-1$ values in the diagonal blocks, and $\tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{N}}}_{rs}=\tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{N}}}_{tu}$ for $r\neq s, t\neq u \in {1,2,...,k}$ and $\sum_i c_i = 0$. Designs that meet this condition include complete randomization, cluster-randomization and block-randomization.
definition[Generalized Neyman variance bound]
The “Generalized Neyman bound" is the is the bound corresponding to the block-diagonal bounding matrix, $\tilde{\mathbf{d}}^{{\scriptscriptstyle{\textnormal{N}}}}$, with block $(r,r)$ given by,
\begin{align*}
\tilde{\mathbf{d}}^{\scriptscriptstyle{N}}_{rr} :=& \sum_{s=1}^{k} \frac{c_r}{c_s}\mathbf{d}_{rs}
\end{align*}
where $c_r$ and $c_s$ are, respectively, elements $r$ and $s$ from from the contrast vector, $c$.
theoremThe generalized Neyman bound, with $\tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{N}}}$ given in Definition (ref) is an identified variance bound when partitions $\textnormal{I}(\tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{N}}}_{rr}==-1)=0_{\scriptscriptstyle n\times n}$, i.e., there are no $-1$ values in the diagonal blocks, $\tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{N}}}_{rs}=\tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{N}}}_{tu}$ for $r\neq s, t\neq u \in {1,2,...,k}$, and $\sum_i c_i = 0$.
proofFirst, with $z_c^{\scriptscriptstyle{\textnormal{HT}}} =\textnormal{\footnotesize diag}(y) \mathds{1} (\mathds{1}' \mathds{1}){}^{-1}$ and letting $y_r$ be the length-$n$ vector of potential outcomes for the $r^{th}$ treatment arm and $\tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{N}}}_{rs}$ be the $r,s$ partition of $\mathbf{d}$, we have
\begin{align*}
n^2 {z^{{\scriptscriptstyle{HT}}}}' \mathbf{d} z^{\scriptscriptstyle{HT}} = & y' diag(c'\mathds{1}) \mathbf{d} diag(c'\mathds{1}) y
\\ = & \sum_{r=1}^k c_r^2 y_r' \mathbf{d}_{rr} y_r +\sum_{r=1}^{k-1} \sum_{s=1}^{k} c_r c_s \left(y_r' \mathbf{d}_{rs} y_s+y_s' \mathbf{d}_{sr} y_r\right)
\end{align*}
Next, define the $r,s$ treatment effect as $\tau_{rs}:=y_r-y_s$ and note that $\mathbf{d}_{12}=\mathbf{d}_{rs}$ for $r\neq s$. Then, by the definition of $\tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{N}}}$,
\begin{align*}
n^2 {z^{{\scriptscriptstyle{HT}}}}' \tilde{\mathbf{d}}^{\scriptscriptstyle{N}} z^{\scriptscriptstyle{\textnormal{HT}}} = & y' \textnormal{ diag}(c'\mathds{1}) \tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{N}}} \textnormal{ diag}(c'\mathds{1}) y
\\ = & \sum_{r=1}^k c_r^2 y_r' \mathbf{d}_{rr} y_r +\sum_{r=1}^{k-1} \sum_{s=r+1}^{k} c_r c_s \left(y_r' \mathbf{d}_{rs} y_r+y_s' \mathbf{d}_{sr} y_s\right)
\\ = & \sum_{r=1}^k c_r^2 y_r' \mathbf{d}_{rr} y_r +\sum_{r=1}^{k-1} \sum_{s=r+1}^{k} c_r c_s \left(y_r' \mathbf{d}_{12} y_r+y_s' \mathbf{d}_{12} y_s\right)
\\ = & \sum_{r=1}^k c_r^2 y_r' \mathbf{d}_{rr} y_r +\sum_{r=1}^{k-1} \sum_{s=r+1}^{k} c_r c_s \left(y_r' \mathbf{d}_{12} (y_s+\tau_{rs})+y_s' \mathbf{d}_{12} (y_r-\tau_{rs}) \right)
\\ = & n^2 {z^{{\scriptscriptstyle{\textnormal{HT}}}}}' \mathbf{d} z^{\scriptscriptstyle{\textnormal{HT}}} + \sum_{r=1}^{k-1} \sum_{s=r+1}^{k} c_r c_s \left(y_r' \mathbf{d}_{12}\tau_{rs}- y_s'\mathbf{d}_{12} \tau_{rs} \right)
\\ = & n^2 {z^{{\scriptscriptstyle{\textnormal{HT}}}}}' \mathbf{d} z^{\scriptscriptstyle{\textnormal{HT}}} + \sum_{r=1}^{k-1} \sum_{s=r+1}^{k} c_r c_s \left(y_r' \mathbf{d}_{12}\tau_{rs}- (y_r-\tau_{rs})'\mathbf{d}_{12} \tau_{rs} \right)
\\ = & n^2 {z^{{\scriptscriptstyle{\textnormal{HT}}}}}' \mathbf{d} z^{\scriptscriptstyle{\textnormal{HT}}} + \sum_{r=1}^{k-1} \sum_{s=r+1}^{k} c_r c_s \tau_{rs}'\mathbf{d}_{12} \tau_{rs}.
\end{align*}
Next, to show that the second term is non-negative, note that $\tau_{rs} =\tau_{rk}-\tau_{sk}$, and write
\begin{align*}
\sum_{r=1}^{k-1} \sum_{s=r+1}^{k} c_r c_s \tau_{rs}'\mathbf{d}_{12} \tau_{rs} = & \frac{1}{2}\sum_{r=1}^{k} \sum_{s=1}^{k} c_r c_s \tau_{rs}'\mathbf{d}_{12} \tau_{rs}
\\ = & \frac{1}{2}\sum_{r=1}^{k} \sum_{s=1}^{k} c_r c_s \left(\tau_{rk}-\tau_{sk}\right)'\mathbf{d}_{12} \left(\tau_{rk}-\tau_{sk}\right) \\ = & \frac{1}{2}\sum_{r=1}^{k} \sum_{s=1}^{k} c_r c_s \left(\tau_{rk}'\mathbf{d}_{12} \tau_{rk}+\tau_{sk}'\mathbf{d}_{12} \tau_{sk} - 2 \tau_{sk}'\mathbf{d}_{12} \tau_{rk}\right)
\\ = & \sum_{r=1}^{k} \sum_{s=1}^{k} c_r c_s \tau_{rk}'\mathbf{d}_{12} \tau_{rk} -\sum_{r=1}^{k} \sum_{s=1}^{k} c_r c_s \tau_{sk}'\mathbf{d}_{12} \tau_{rk}
\\ = & \sum_{r=1}^{k} c_r \tau_{rk}'\mathbf{d}_{12} \tau_{rk} \left(\sum_{s=1}^{k} c_s\right) - \left(\sum_{s=1}^{k} c_s \tau_{sk}\right)'\mathbf{d}_{12} \left(\sum_{r=1}^{k}c_r \tau_{rk} \right)
\\ =& 0 - {\tau^*} '\mathbf{d}_{12} \tau^*
\\ \geq& 0
\end{align*}
where the second to last line uses $\sum_{s=1}^{k} c_s=0$ and the definition ${\tau^*}:=\sum_{s=1}^{k} c_s \tau_{sk}$. The last line follows because $\mathbf{d}_{12}$ is negative semidefinite.
A novel proof of the Aronow-Samii bound
Consider an identified bound proposed by aronowsamii17 that has the a unusual virtue of being perfectly general, i.e., applicable to arbitrary (identified) designs.
definition[Aronow-Samii variance bound]
The “Aronow-Samii variance bound" is the bound corresponding to the bounding matrix,
\begin{align*}
\tilde{\mathbf{d}}^{\scriptscriptstyle{AS}}:=\mathbf{d} +I\left(\mathbf{d}=-1\right)+ diag(
I\left(\mathbf{d}=-1\right) 1_{\scriptscriptstyle {kn}}
)
\end{align*}
where the indicator function, $\textnormal{I}(\mathbf{d}=-1)$, returns a matrix of with ones indicating the location of -1 entries in $\mathbf{d}$ and zeros elsewhere, and $\textnormal{\footnotesize diag}(.)$ creates a diagonal matrix from a vector.
theoremThe Aronow-Samii variance bound, $ n^{-2} y' \tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{AS}}} y$, is an identified bound for $ n^{-2} y' {\mathbf{d}} y$.
proofBy definition of $\tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{AS}}}$,
\begin{align*}
\tilde{\mathbf{d}}^{\scriptscriptstyle{AS}}-\mathbf{d}= I\left(\mathbf{d}=-1\right)+\scriptsize diag\left(I\left(\mathbf{d}=-1\right) 1_{\scriptscriptstyle kn}\right).
\end{align*}
Note that by construction ($\tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{AS}}}-\mathbf{d}$) has diagonal elements set equal to the sum of the off-diagonal elements in its row (which by construction are either 0 or 1). The Gershgorin circle theorem implies that a real matrix is positive semi-definite if, for all $i$, the $i^{th}$ diagonal element is greater or equal to the sum of the absolute values of the other elements in the $i^{th}$ row. So, by the Gershgorin circle theorem $\tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{AS}}}-\mathbf{d}$ is positive semidefinite. Therefore, by Lemma ((ref)), $n^{-2} y'\tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{AS}}} y$ is a variance bound. Moreover, as long as the design is an identified design (i.e., $0<\pi_{1i}<1$ for all $i$), it is an identified bound because $\text{\textnormal{I}}\left(\mathbf{d}=-1\right)$ ensures that the elements of $\mathbf{d}$ equal to $-1$ correspond to 0's in $\tilde{\mathbf{d}}^{{\scriptscriptstyle{\textnormal{AS}}}}$.
remarkAronow and Samii (2017) derive their bound using Young's inequality. The above-theorem and proof using the Gershgorin circle theorem tie their insight to the current framework.
Proposed algorithm for variance bounds for any design
The following is an algorithm which that can obtain an identified variance bound. Like the AS bound it has the virtue of being applicable to virtually any design. The algorithm is a proof of concept, demonstrating the utility of the notation scheme which allows for the application of matrix theory for the creation of alternative bounds. The subject of comparing bounds will be considered further in Section (ref).
algorithm[algorithm omitted — 33 chars of source]
enumerate• Initialize $kn \times kn$ matrix $\mathbf{t}$. Examples could be $\text{\textnormal{I}}(\mathbf{d}=-1)$ or, if the conditions for the Neyman bound not be applicable, start with $\tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{N}}}-\mathbf{d}$ which may approximate a bound
• Obtain the eigen decomposition of matrix $\mathbf{t}$. If all eigenvalues are non-negative (within tolerance), goto Step 6, otherwise continue
• Update $\mathbf{t}=\mathbf{v} (\mathbf{e} \circ \text{\textnormal{I}}(\mathbf{e}>0)) \mathbf{v}'$ where $\mathbf{v}$ is the matrix of eigenvectors and $\mathbf{e}$ is a diagonal matrix of eigenvalues
• Update $\mathbf{t}=\text{\textnormal{I}}(\mathbf{d}=-1)+\text{\textnormal{I}}(\mathbf{d}\neq -1) \circ \mathbf{t}$
• Return to Step 2
• Set $\tilde{ \mathbf{d} }^{\scriptscriptstyle{\textnormal{M}}} =\mathbf{d}+\mathbf{t}$
As above, $\circ$ is elementwise multiplication and, for example, $\text{\textnormal{I}}(\mathbf{e}>0)$ is an indicator function returning a matrix of ones and zeros indicating which elements of $\mathbf{e}$ are greater than zero.
Conceptually, the goal of the algorithm is to create a matrix $\mathbf{t}$ that can be added to $\mathbf{d}$ yielding a $\tilde{ \mathbf{d} }$ matrix that corresponds to an identified variance bound. By Lemma (ref) and Definition (ref), there are two requirements for $\mathbf{t}$. First it must be positive semi-definite, and, second, elements corresponding to $-1$'s in the matrix $\mathbf{d}$ must equal one. In step 1, $\mathbf{t}$ meets the second criterion, but not the first. In step 3, the algorithm creates an approximation to the initial $\mathbf{t}$ matrix by way of the eigen decomposition that ensures positive semi-definiteness, thus meeting the first criterion. However, due to the approximation, $\mathbf{t}$ no longer meets the second criterion. Therefore, in step 4 the algorithm forces $\mathbf{t}$ to have 1's wherever $\mathbf{d}$ has $-1$'s in order to again meet the second criteria. But doing so means that $\mathbf{t}$ will no longer meet the first criteria. So, the algorithm iterates through steps 2-4 until convergence is achieved (i.e., until all eigenvalues are non-negative in step 2) at which point $\mathbf{t}$ meets both criteria and, thus, $\tilde{ \mathbf{d} }^{\scriptscriptstyle{\textnormal{M}}}$ corresponds to an identified bound.
Comparing bounds
definition[Tighter bound]
Let $\tilde{\mathbf{d}}^a$ and $\tilde{\mathbf{d}}^b$ correspond to two identified bounds. Matrix $\tilde{\mathbf{d}}^a$ is corresponds to a tighter bound than $\tilde{\mathbf{d}}^b$ if $\tilde{ \mathbf{d} }^b - \tilde{ \mathbf{d} }^a$ is positive semidefinite.
definition[Invariant bounding matrix]
A matrix $\tilde{\mathbf{d}}$ is an invariant bounding matrix if it is an bounding matrix and if all $n \times n$ partitions, $\tilde{\mathbf{d}}_{ij} 1_{\scriptscriptstyle {kn}} = 0_{\scriptscriptstyle kn}$, i.e., all rows of the partition (or, equivalently, all columns) sum to zero.
Example (Paired randomization) : Consider a pair-randomized design, whereby units are “blocked" (i.e., stratified) into groups of two, and then, in each block, one of the two units is randomly assigned to treatment while the other is assigned to control. Assignments across blocks are independent.
When $n=4$ (and assuming w.l.o.g. that the data are sorted by pair), the matrix $\mathbf{d}$ is
align*[align* omitted — 506 chars of source]
noting that empty cells represent 0.
For the pair-randomized design, the Neyman bound cannot be applied because $\mathbf{d}_{00}$ and $\mathbf{d}_{11}$ have negative entries. The Aronow-Samii bound and Algorithm (ref) have bounding matrices
center[center omitted — 1,046 chars of source]
respectively. By the Gershgorian circle theorem the difference,
align*[align* omitted — 465 chars of source]
is positive semi-definite, proving that $\tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{M}}}$ corresponds to a tighter variance bound. Confirmation also comes from eigendecomposition of the difference, $\tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{AS}}} - \tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{M}}}$, which yields all non-negative eigenvalues: 2, 2, 2, 2, 0, 0, 0, and 0.
One might alternatively choose the invariant bounding matrix,
align*[align* omitted — 605 chars of source]
The bound can be verified because the eigenvalues of $\tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{INVAR}}}-\mathbf{d}$ are 8, 0, 0, 0, 0, 0, 0, and 0. However, eigendecomposition of $\tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{INVAR}}}-\tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{M}}}$ gives eigenvalues 4, 0, 0, 0, 0, 0, 0, and -4, indicating that the better bound may depend on outcome vector, $y$, and perhaps the estimator as well.
{ $\triangle$}
Variance bound estimation
With an identified variance bounds defined and several methods of obtaining matrices, $\tilde{\mathbf{d}}$, this section turns to the subject of variance bound estimation.
First define the $kn \times kn$ matrix of probabilities and joint probabilities of assignment,
align[align omitted — 150 chars of source]
Next define an inverse probability weighted version of bounding matrix, $\tilde{ \mathbf{d} }$, as
align[align omitted — 141 chars of source]
with $/$ denoting element-wise division defined such that division by zero equals zero. Then an unbiased estimator of a variance bound for the Horvitz-Thompson estimator can be written,
align[align omitted — 345 chars of source]
with $z_c^{\scriptscriptstyle{\textnormal{HT}}} := \textnormal{\footnotesize diag}(y)\mathds{1} \left(\mathds{1}' \mathds{1} \right)^{-1} c$.
It is unbiased for the variance bound ${z_c^{\scriptscriptstyle{\textnormal{HT}}}}' \tilde{\mathbf{d}} z_c^{\scriptscriptstyle{\textnormal{HT}}}$ because $\text{\textnormal{E}} \left[\mathbf{R} \tilde{ \mathbf{d} }_{\hspace{-.6mm}{}{/}}{}_{\scriptscriptstyle \hspace{-.6mm}\mathbf{p}} \mathbf{R} \right]=\tilde{ \mathbf{d} }$ by construction. Being inverse-probability weighted, the variance bound estimator in ((ref)) is, itself, a Horvitz-Thompson estimator.
For other linear estimators, examples of which are given in Table (ref), the bound $z_c' \tilde{ \mathbf{d} }z_c$ cannot be estimated unbiasedly because the definition of $z_c$ will often include quantities that, themselves, must be estimated. However, an appeal to the plug-in principle suggests the use of
align[align omitted — 288 chars of source]
with $\widehat{z}_c$ having the same form as $z_c$ but with sample analogues replacing some components.
Example (The special case of Eicker-Huber-White (a.k.a. “heteroskedastic consistent", “sandwich", and “robust") standard errors) : For the OLS estimator, $z^{\scriptscriptstyle{\textnormal{OLS}}}_c$ is defined in Table ((ref)). The plug-in principle motivates the use of
align*[align* omitted — 233 chars of source]
where $\mathbf{R} \widehat{u}:=\mathbf{R} (y-\mathbb{x} \widehat{b}^{{\scriptscriptstyle{\textnormal{OLS}}}} )$ and $\widehat{b}^{{\scriptscriptstyle{\textnormal{OLS}}}}:=\left( \mathbb{x}' \mathbf{R} \mathbb{x}\right)^{-1} \mathbb{x}' \mathbf{R} y$ is the OLS coefficient. Then from equation ((ref)) we have,
align*[align* omitted — 582 chars of source]
This is a the variance bound estimator in ((ref)) made specific to OLS. So far it is applicable to virtually any design and any variance bound.
Next, specify a Bernoulli design, in which units are assigned independently to treatment. (Probabilities of assignment may be equal across units, but they need not be in this example.) In this design, the diagonal elements of $\mathbf{d}$ are equal to the diagonal of $\boldsymbol{\pi}^{-1} -\mathbf{i}_{\scriptscriptstyle kn}$, where $\mathbf{i}_{\scriptscriptstyle kn}$ is an identity matrix. Further, any of the above bounding methods yields $\tilde{\mathbf{d}}=\boldsymbol{\pi}^{-1} -\mathbf{i} + \mathbf{i}=\boldsymbol{\pi}^{-1}$. Thus $\tilde{ \mathbf{d} }_{\hspace{-.6mm}{}{/}}{}_{\scriptscriptstyle \hspace{-.6mm}\mathbf{p}}=\boldsymbol{\pi}^{-2}$ so that $\boldsymbol{\pi} \tilde{ \mathbf{d} }_{\hspace{-.6mm}{}{/}}{}_{\scriptscriptstyle \hspace{-.6mm}\mathbf{p}} \boldsymbol{\pi}=\mathbf{i}_{\scriptscriptstyle kn}$ is the identity matrix. So the OLS variance bound estimator for Bernoulli designs simplifies to,
align*[align* omitted — 414 chars of source]
This is White's (1980) canonical “sandwich" variance estimator, sometimes referred to as HC0. { $\triangle$}
remarkThe example shows that Eicker-Huber-White standard errors are a special case of ((ref)) for OLS in a Bernoulli design. Note, however, that ((ref)) is much more general. It applies to any linear estimator, virtually any design and any (identified) variance bound.
remarkAdjustments for degrees of freedom (e.g., HC1) or leverage (e.g., HC2, HC3, etc.) can be applied as well.
Example (The special case of “cluster robust" standard errors) : Also consider this variance bound estimator for OLS in designs in which clusters are assigned independently to treatment. Then, if we choose the Neyman bound $\tilde{ \mathbf{d} }_{\hspace{-.6mm}{}{/}}{}_{\scriptscriptstyle \hspace{-.6mm}\mathbf{p}}^{\scriptscriptstyle{\textnormal{N}}}$ (or $\tilde{ \mathbf{d} }_{\hspace{-.6mm}{}{/}}{}_{\scriptscriptstyle \hspace{-.6mm}\mathbf{p}}^{\scriptscriptstyle{\textnormal{M}}}$, which is equivalent in the case of for Bernoulli assignment of clusters), and assuming w.l.o.g. that units are sorted by cluster, then $\boldsymbol{\pi} \tilde{ \mathbf{d} }_{\hspace{-.6mm}{}{/}}{}_{\scriptscriptstyle \hspace{-.6mm}\mathbf{p}}^{{\scriptscriptstyle{\textnormal{N}}}} \boldsymbol{\pi} $ resolves to a block diagonal matrix of 1's with the blocks corresponding to clusters. Hence, ((ref)) also reproduces the “cluster-robust" standard errors sometimes referred to as CR0 as a special case.{ $\triangle$}
Asymptotics
Conditions for Convergence of Horvitz-Thompson Estimators
First establishing the unbiasedness of Horvitz-Thompson estimators will allow for straightforward proofs of consistency.
lemmaThe Horvitz-Thompson estimator for an outcome vector, $y$, and given contrast, $c$, is unbiased for $\delta_c$.
proof\begin{align*}
E \left[ c' \left(\mathds{1}' \mathds{1} \right)^{-1} \mathds{1}' \boldsymbol{\pi}^{-1} \mathbf{R} y \right]= & c' \left(\mathds{1}' \mathds{1} \right)^{-1} \mathds{1}' \boldsymbol{\pi}^{-1} E \left[ \mathbf{R} \right]y
\\ = & c' \left(\mathds{1}' \mathds{1} \right)^{-1} \mathds{1}' \boldsymbol{\pi}^{-1} \boldsymbol{\pi} y
\\ = &c' \left(\mathds{1}' \mathds{1} \right)^{-1} \mathds{1}'y
\\ = & \delta_c,
\end{align*}
condition[Bounded contrast]
The chosen contrast vector is finite, i.e., there exists a finite value $u_c$ such that $\textnormal{max}(|c|)< u_c$.
condition[Bounded outcomes]
There exists a finite value, $u_y$, such that $\textnormal{max}(|y|)<u_y$ for all $n$.
condition[Design constraint for consistent Horvitz-Thompson estimators]
There exists a finite value, $u_\mathbf{d}$, such that $n^{-1} ||\mathbf{d}||_{{1,1}} < u_\mathbf{d}$ for all $n$, where $||.||_{1,1}$ is the matrix norm that sums the absolute values of the matrix entries.
theorem[Root-n consistency of HT estimators]
By Lemma (ref) and Conditions (ref)-(ref) the Horvitz-Thompson estimator is root-n consistent.
proofGiven Lemma (ref), it is sufficient to show that the variance converges at the parametric rate, i.e., that $n \text{\textnormal{V}}(\delta_c^{\scriptscriptstyle{\textnormal{HT}}})$ is bounded, in order to prove consistency. By Holder's Inequality,
\begin{align*}
n V(\delta_c^{\scriptscriptstyle{HT}}) & \leq n max(|c'\left(\mathds{1}' \mathds{1} \right)^{-1}\mathds{1}'y|)^2 ||\mathbf{d}||_{{1,1}}
\\ & \leq max(|c|)^2 max(|y|)^2 n^{-1} ||\mathbf{d}||_{{1,1}}
\\ & \leq u_c^2 u_y^2 u_\mathbf{d} ,
\end{align*}
with the last line using Conditions (ref)-(ref).
Example (Checking consistency of HT estimators for completely randomized experiments) : Consider a completely randomized experiment with $n$ units where a fixed number of units, $n_c$, are randomly assigned to control, and the remainder, $n_t= n - n_c$, are assigned to treatment. Assume an asymptotic sequence of designs is such that there exists a constant value, $\pi$, such that $\frac{n_t}{n} \rightarrow \pi$ as $n \rightarrow \infty$ with $0<\pi_t<1$. Partition $\mathbf{d}$ into four $(n \times n)$ matrices and let $\mathbf{d}_{ab}$ represent the $a,b \in \{1,2\}$ partition. Each partition has elements which take on two possible values, one on the diagonal and another on the off-diagonal. Because $\mathbf{d}_{12}=\mathbf{d}_{21}$, entries of matrix $\mathbf{d}$ take on one of six possible values. In Table (ref), analysis of the these six values and their corresponding frequencies shows that a completely randomized design yields $\frac{1}{n}||\mathbf{d}||_{{1,1}}=2\left(\frac{n_t}{n_c}+\frac{n_c}{n_t}+2\right)=O(1)$. Thus, Condition (ref) is satisfied. Therefore, by Theorem (ref), Horvitz-Thompson estimators are consistent for completely randomized experiments for bounded contrast, $c$, and outcome vector, $y$. { $\triangle$}
table[table omitted — 1,298 chars of source]
Conditions for Convergence of WLS estimator class
condition[Bounded covariates]
There exists a finite value $u_\mathbf{x}$ that bounds the covariate values, i.e., $\textnormal{max}(|\mathbf{x} |)<u_\mathbf{x} $, for all $n$.
condition[Bounded $\boldsymbol{\pi} \mathbf{m}$]
There exists a finite value $u_{\boldsymbol{\pi}\mathbf{m}}$ that bounds $\boldsymbol{\pi}$ times the WLS “weighting" matrix $\mathbf{m}$, i.e., $\textnormal{max}(| \boldsymbol{\pi} \mathbf{m} |)<u_{\boldsymbol{\pi}\mathbf{m}} $, for all $n$.
lemma[Root-n consistency of WLS]
By conditions (ref)-(ref) and Theorem (ref), the WLS “numerator” vector, $\frac{1}{n} \mathbb{x}' \mathbf{m} \mathbf{R} y$, is root-n consistent for $\frac{1}{n} \mathbb{x}' \mathbf{m} \boldsymbol{\pi} y$. Likewise, by conditions (ref)-(ref) and Theorem (ref), the WLS “denominator” matrix, $\frac{1}{n} \mathbb{x}' \mathbf{m} \mathbf{R} \mathbb{x}$, is root-n consistent for $\frac{1}{n} \mathbb{x}' \mathbf{m} \boldsymbol{\pi} \mathbb{x}$. Further, by the continuous mapping theorem $\widehat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}} \rightarrow b^{\scriptscriptstyle{\textnormal{WLS}}}$.
proofLet $\mathbb{x}_i$ be the column vector created from the $i^{th}$ column of $\mathbb{x}$. Then the $i^{the}$ element of $\frac{1}{n} \mathbb{x} \mathbf{m} \mathbf{R} y $ can be written,
\begin{align*}
\{ \frac{1}{n} \mathbb{x} \mathbf{m} \mathbf{R} y \}_{i } =& \frac{1}{n} \mathbb{x}_i' \mathbf{m} \mathbf{R} y
\\ = & \frac{1}{n} 1_{\scriptscriptstyle {kn}}' diag\left(\mathbb{x}_i\right) \mathbf{m} \mathbf{R} y
\\ = & \frac{1}{n} 1_{\scriptscriptstyle {kn}}' \mathbf{R} \mathbf{m} diag\left(\mathbb{x}_i\right) y
\\ = & 1_{\scriptscriptstyle {k}}' \mathbf{w}^{{\scriptscriptstyle{HT}}} \mathbf{R} \boldsymbol{\pi} \mathbf{m} diag\left(\mathbb{x}_i\right) y
\\ = & 1_{\scriptscriptstyle {k}}' \mathbf{w}^{{\scriptscriptstyle{HT}}} \mathbf{R} q
\end{align*}
with $q=\boldsymbol{\pi} \mathbf{m} \textnormal{\footnotesize diag}\left(\mathbb{x}_i\right) y $, showing that the elements of the denominator matrix are Horvitz-Thompson estimators with outcome vector $q$ and contrast vector $ 1_{\scriptscriptstyle {k}}$. Now by Theorem (ref), this is consistent because $q$ is bounded, i.e., $\textnormal{max}(|q|) \leq u_{\boldsymbol{\pi} \mathbf{m}} u_\mathbf{x} u_y $.
Similarly, the $i,j$ element of the WLS “denominator” matrix, can be written
\begin{align*}
\{ \frac{1}{n} \mathbb{x} \mathbf{m} \mathbf{R} \mathbb{x} \}_{ij} = & 1_{\scriptscriptstyle {k}}' \mathbf{w}^{{\scriptscriptstyle{HT}}} \mathbf{R} r
\end{align*}
with $r=\boldsymbol{\pi} \mathbf{m} \textnormal{\footnotesize diag}\left(\mathbb{x}_i\right) \mathbb{x}_j $, showing that the elements of the denominator matrix are Horvitz-Thompson estimators with outcome vector $r$ and contrast vector $ 1_{\scriptscriptstyle {k}}$. Now by Theorem (ref), this is consistent because $r$ is bounded, i.e., $\textnormal{max}(|r|) \leq u_{\boldsymbol{\pi} \mathbf{m}} u_\mathbf{x}^2 $.
condition[Stability of WLS “denominator" estimand]
The denominator of the “true" WLS coefficient, $\frac{1}{n} \mathbb{x}' \mathbf{m} \pi \mathbb{x}$, is invertable for all $n$ and converges in probability to a matrix, $\mathbf{v}$, with finite entries.
theorem[Consistency of the Taylor approximation]
By Conditions (ref)-(ref), the Taylor approximate coefficient, $ \widehat{b}^{{\scriptscriptstyle{\textnormal{T}}}({\scriptscriptstyle{\textnormal{WLS}}})}:= b^{{\scriptscriptstyle{\textnormal{WLS}}}} + \left(\mathbb{x}' \mathbf{m} \boldsymbol{\pi} \mathbb{x} \right)^{-1}\mathbb{x}' \mathbf{m} \mathbf{R} (y-\mathbb{x} b^{{\scriptscriptstyle{\textnormal{WLS}}}})$ is root-n consistent for ${b}^{{\scriptscriptstyle{\textnormal{WLS}}}}:= \left(\mathbb{x}' \mathbf{m} \boldsymbol{\pi} \mathbb{x} \right)^{-1}\mathbb{x}' \mathbf{m} \boldsymbol{\pi} y$, i.e., $ \widehat{b}^{{\scriptscriptstyle{\textnormal{T}}}({\scriptscriptstyle{\textnormal{WLS}}})}- b^{{\scriptscriptstyle{\textnormal{WLS}}}} = O_p(1/\sqrt{n})$.
proofWe have
\begin{align*}
\widehat{b}^{{\scriptscriptstyle{T}}({\scriptscriptstyle{WLS}})}- b^{{\scriptscriptstyle{WLS}}}:= &\left(\mathbb{x}' \mathbf{m} \boldsymbol{\pi} \mathbb{x} \right)^{-1}\mathbb{x}' \mathbf{m} \mathbf{R} (y-\mathbb{x} b^{{\scriptscriptstyle{WLS}}})
\\ =& \left(\mathbb{x}' \mathbf{m} \boldsymbol{\pi} \mathbb{x} \right)^{-1} \left( \left(\mathbb{x}' \mathbf{m}\mathbf{R} y-\mathbb{x}' \mathbf{m} \boldsymbol{\pi} y\right) -\left( \mathbb{x}' \mathbf{m} \mathbf{R} \mathbb{x} -\mathbb{x}' \mathbf{m} \boldsymbol{\pi} \mathbb{x} \right) b^{{\scriptscriptstyle{WLS}}} \right)
\\ =& \left(\frac{1}{n} \mathbb{x}' \mathbf{m} \boldsymbol{\pi} \mathbb{x} \right)^{-1} \left( \left(\frac{1}{n} \mathbb{x}' \mathbf{m}\mathbf{R} y-\frac{1}{n}\mathbb{x}' \mathbf{m} \boldsymbol{\pi} y\right) -
\left( \frac{1}{n}\mathbb{x}' \mathbf{m} \mathbf{R} \mathbb{x} - \frac{1}{n} \mathbb{x}' \mathbf{m} \boldsymbol{\pi} \mathbb{x} \right) b^{{\scriptscriptstyle{WLS}}} \right)
\\ =& O_p(1) \left( O_p(1/ \sqrt{n})+O_p(1/ \sqrt{n}) \right),
\\=& O_p(1/ \sqrt{n})
\end{align*}
where the second to last line uses Lemma (ref) and Condition (ref).
theorem[Asymptotic Equivalence of WLS and its Taylor Approximation]
By Lemma (ref) and Condition (ref), WLS is asymptotically equivalent to the Taylor linear approximation for WLS.
proofLet the “true" WLS coefficient be $b^{\scriptscriptstyle{\textnormal{WLS}}} = \left(\mathbb{x}' \mathbf{m} \boldsymbol{\pi} \mathbb{x} \right)^{-1}\mathbb{x}' \mathbf{m} \boldsymbol{\pi} y$, then write the WLS estimator as,
\begin{align*}
c' \mathbf{W}^{\scriptscriptstyle{WLS}} \mathbf{R} y = & c' \left( \mathbb{x}' \mathbf{m} \mathbf{R} \mathbb{x} \right)^{-1} \mathbb{x}' \mathbf{m} \mathbf{R} y
\\ = & c' b^{{\scriptscriptstyle{WLS}}} + c' \left( \mathbb{x}' \mathbf{m} \mathbf{R} \mathbb{x} \right)^{-1} \mathbb{x}' \mathbf{m} \mathbf{R} \left( y - \mathbb{x} b^{{\scriptscriptstyle{WLS}}}\right)
\\ = & c' b^{{\scriptscriptstyle{WLS}}} + c' \left( \mathbb{x}' \mathbf{m} \boldsymbol{\pi} \mathbb{x} \right)^{-1} \mathbb{x}' \mathbf{m} \mathbf{R} \left( y - \mathbb{x} b^{{\scriptscriptstyle{WLS}}}\right)
\\ & + c' \left( \left( \mathbb{x}' \mathbf{m} \mathbf{R} \mathbb{x} \right)^{-1} -\left( \mathbb{x}' \mathbf{m} \boldsymbol{\pi} \mathbb{x} \right)^{-1} \right)\mathbb{x}' \mathbf{m} \mathbf{R} \left( y - \mathbb{x} b^{{\scriptscriptstyle{WLS}}}\right)
\\ = & c' b^{{\scriptscriptstyle{\textnormal{WLS}}}} + c' \left( \mathbb{x}' \mathbf{m} \boldsymbol{\pi} \mathbb{x} \right)^{-1} \mathbb{x}' \mathbf{m} \mathbf{R} \left( y - \mathbb{x} b^{{\scriptscriptstyle{\textnormal{WLS}}}}\right)
\\ & + c' \left( \left( \mathbb{x}' \mathbf{m} \mathbf{R} \mathbb{x} \right)^{-1} -\left( \mathbb{x}' \mathbf{m} \boldsymbol{\pi} \mathbb{x} \right)^{-1} \right)\mathbb{x}' \mathbf{m} \mathbf{R} \left( \mathbb{x} \widehat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}} - \mathbb{x} b^{{\scriptscriptstyle{\textnormal{WLS}}}}\right)
\\ & + c' \left( \left( \mathbb{x}' \mathbf{m} \mathbf{R} \mathbb{x} \right)^{-1} -\left( \mathbb{x}' \mathbf{m} \boldsymbol{\pi} \mathbb{x} \right)^{-1} \right)\mathbb{x}' \mathbf{m} \mathbf{R} \left( y- \mathbb{x} \widehat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}} \right)
\\ = & \widehat{\delta}_c^{{\scriptscriptstyle{\textnormal{T}}}({\scriptscriptstyle{\textnormal{WLS}}})} + c' \left( \left( \frac{1}{n}\mathbb{x}' \mathbf{m} \mathbf{R} \mathbb{x} \right)^{-1} -\left( \frac{1}{n} \mathbb{x}' \mathbf{m} \boldsymbol{\pi} \mathbb{x} \right)^{-1} \right)\left(\frac{1}{n}\mathbb{x}' \mathbf{m} \mathbf{R} \mathbb{x} \right) \left( \widehat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}} - b^{{\scriptscriptstyle{\textnormal{WLS}}}}\right)
\\ = & \widehat{\delta}_c^{{\scriptscriptstyle{\textnormal{T}}}({\scriptscriptstyle{\textnormal{WLS}}})} + O_p(1/\sqrt{n})O_p(1)O_p(1/\sqrt{n})
\\ = & \widehat{\delta}_c^{{\scriptscriptstyle{\textnormal{T}}}({\scriptscriptstyle{\textnormal{WLS}}})} + O_p(1/n)
\end{align*}
Conditions for consistent variance estimation
definition[Second-order design matrix]
The “second-order design matrix" is a forth order tensor $(kn \times kn \times kn \times kn)$ of variances and covariances of inverse-probability weighted pairwise-joint inclusion indicators, written,
\begin{align*}
\mathbf{d}{ l} :=\Big( E\left[\left(\mathbf{R} 1_{\scriptscriptstyle n}1_{\scriptscriptstyle kn}'\mathbf{R}\right) \otimes \left( \mathbf{R} 1_{\scriptscriptstyle kn}1_{\scriptscriptstyle kn}'\mathbf{R}\right) \right] -\mathbf{p} \otimes \mathbf{p} \Big) / \left( \mathbf{p} \otimes \mathbf{p} \right),
\end{align*}
where $\mathbf{p}:=\text{\textnormal{E}}\left[\mathbf{R} 1_{\scriptscriptstyle n}1_{\scriptscriptstyle kn}'\mathbf{R}\right]$ is a matrix of with inclusion probabilities on the diagonal and pair-wise joint inclusion probabilities off the diagonal, “$\otimes$" is the tensor outer product and “/" is elementwise division with division by zero resolving to zero.
condition[Second order design constraint for consistent variance estimation]
There exists a finite constant $u_{ \mathbf{d}{ \hspace{-.6 mm} \textnormal{l}}}$ such that $\frac{1}{n} \left|\left|\left(\tilde{\mathbf{d}} \otimes \tilde{\mathbf{d}} \right) \circ \mathbf{d}{ \hspace{-.6 mm} \textnormal{l}} \right|\right|_{1,1,1,1} < u_{ \mathbf{d}{ \hspace{-.6 mm} \textnormal{l}}}$ for all $n$,
where “$\otimes$" is tensor outer product, “$\circ$" is elementwise multiplication, $\left|\left|. \right|\right|_{1,1,1,1}$ gives the sum of the absolute values of the tensor entries.
theorem[Consistency of the Horvitz-Thompson variance estimator]
By Conditions (ref)-(ref) and (ref) the variance estimator for the Horvitz-Thompson point estimator is consistent.
proofThe variance of $n$ times the Horvitz-Thompson variance estimator (times $n$) is,
\begin{align*}
n V \bigg ( n \widehat{\tilde{V}}\left(\widehat{\delta}^{{\scriptscriptstyle{HT}}} \right) \bigg ) = & n\textnormal{E} \left[ \left(n {z^{\scriptscriptstyle{\textnormal{HT}}}}' \mathbf{R} \tilde{ \mathbf{d} }_{{/}}_{\scriptscriptstyle \mathbf{p}} \mathbf{R} z^{\scriptscriptstyle{\textnormal{HT}}} - n {z^{\scriptscriptstyle{\textnormal{HT}}}}' \tilde{\mathbf{d}} z^{\scriptscriptstyle{\textnormal{HT}}} \right)^2\right]
\\ \leq & \max\left(\left|c\right|\right)^4 \max\left(\left|y\right|\right)^4 \frac{1}{n} \left| \left|\left(\tilde{\mathbf{d}} \otimes \tilde{\mathbf{d}} \right) \circ \mathbf{d}{ \textnormal{l}} \right|\right|_{1,1,1,1}
\\ \leq & u_c^4 u_y^4 u_{ \mathbf{d}{ \textnormal{l}}},
\end{align*}
where the second line uses Holder's inequality and the last line uses Conditions (ref), (ref) and (ref).
Example (Checking consistency of HT variance (bound) estimator for completely randomized experiments) : Again consider a completely randomized experiment with $n$ units where a fixed number of units, $n_c$, are randomly assigned to control, and the remainder, $n_t= n - n_c$, are assigned to treatment. Assume an asymptotic sequence of designs is such that there exists a constant value, $\pi$, such that $\frac{n_t}{n} \rightarrow \pi$ as $n \rightarrow \infty$ with $0<\pi<1$. Let the variance bound be the Neyman bound, i.e., $\tilde{\mathbf{d}} =\tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{N}}}$ from Definition (ref). Table (ref) enumerates the unique values that appear in $\frac{1}{n} \left|\left|\left(\tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{N}}} \otimes \tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{N}}} \right) \circ \mathbf{d}{ \hspace{-.6 mm} \textnormal{l}} \right|\right|_{1,1,1,1}$ along with their relative frequencies and shows that Condition (ref) is satisfied for completely randomized experiments. Hence, the Horvitz-Thompson variance (bound) estimator given in Equation (ref) is consistent for the Neyman bound. { $\triangle$}
landscape\begin{table}
\begin{tabular}{|c|c|c|c|c|c| c|} \hline
\makecell{Treatment (T) \\ or Control (C)} & $ijkl$ pattern & count & $\makecell{ \hspace{-8mm}\{ \scriptstyle \text{\textnormal{E}}\left[\left(\mathbf{R} 1_{\scriptscriptstyle n}1_{\scriptscriptstyle kn}'\mathbf{R}\right) \right. \\ \left. \scriptstyle \otimes \left(\mathbf{R} 1_{\scriptscriptstyle n}1_{\scriptscriptstyle kn}'\mathbf{R}\right)\right] \}_{\scriptscriptstyle ijkl}}$ & $\{ \mathbf{p} \otimes \mathbf{p}\}_{\scriptscriptstyle ijkl}$ & $\{\tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{N}}} \otimes \tilde{\mathbf{d}}^{\scriptscriptstyle{\textnormal{N}}} \}_{\scriptscriptstyle ijkl}$
& count $\times \scriptstyle \frac{1}{n} \big \{ \left(\tilde{\mathbf{d}} \otimes \tilde{\mathbf{d}} \right) \circ \mathbf{d}{ \hspace{-.6 mm} \textnormal{l}} \big \} _{\scriptscriptstyle ijkl}
$
\\ \hline $\scriptstyle i,j,k,l \in \textnormal{C}$ & $\scriptstyle i=j=k=l$ & $\scriptstyle n$
& $\frac{n_c}{n}$ & $\frac{n_c^2}{n^2}$
& $\frac{n^2}{n_c^2}$ & $\scriptstyle \frac{n^2 n_t }{n_c^3}=O(1)$
\\ & \makecell{$\scriptstyle i=j=k,l$ or $\scriptstyle i=j=l,k$ \\ or $\scriptstyle i,j=k=l$ or $\scriptstyle i=k=l,j$} & $\scriptstyle 4n(n-1)$ & $\frac{n_c (n_c-1)}{n(n-1)}$ & $\frac{n_c^2(n_c-1)}{n^2(n-1)}$
& $-\frac{n^2}{n_c^2(n-1)}$ & $\scriptstyle -\frac{4n^2 n_t}{n_c^3}=O(1)$
\\ & $\scriptstyle i=j, k=l$ & $\scriptstyle n(n-1)$ & $ \frac{n_c (n_c-1)}{n(n-1)}$ & $\frac{n_c^2 }{n^2 }$
& $\frac{n^2}{n_c^2}$ & $\scriptstyle \frac{n_t n^2}{n_c^3}=O(1)$
\\ & $\scriptstyle i=j,k,l$ or $\scriptstyle i,j,k=l$ & $\scriptstyle 2n(n-1)(n-2)$ & $ \frac{n_c (n_c-1)(n_c-2)}{n(n-1)(n-2)}$ & $\frac{n_c^2(n_c-1) }{n^2(n-1) }$
& $-\frac{n^2}{n_c^2(n-1)}$ & $\scriptstyle \frac{4n_t n^2}{ n_c^3}=O(1)$
\\ & $\scriptstyle i=k,j=l$ or $\scriptstyle i=l, j=k$ & $\scriptstyle 2n(n-1)$ &
$ \frac{n_c (n_c-1)}{n(n-1)}$ & $\frac{n_c^2(n_c-1)^2 }{n^2(n-1)^2 }$
& $-\frac{n^2}{n_c^2(n-1)^2}$ & $\scriptstyle -\frac{2 n^2}{n_c^2}=O(1)$
\\ & \makecell{$\scriptstyle i=k, j, l$ or $\scriptstyle i=l,j,k$ \\ or $\scriptstyle i,j=k,l$ or $\scriptstyle i,j=l,k$} & $\scriptstyle 4n (n-1)(n-2)$
& $ \frac{n_c (n_c-1)(n_c-2)}{n(n-1)(n-2)}$ & $\frac{n_c^2(n_c-1)^2 }{n^2(n-1)^2 }$
& $\frac{n^2}{n_c^2(n-1)^2}$ & $\scriptstyle \frac{4 n^2}{n_c^2}\left(\frac{n(n_c-2)}{n_c(n_c-1)}-\frac{(n-2)}{(n-1)} \right)=O(1)$
\\ & $\scriptstyle i,j,k,l$ & $\scriptscriptstyle n(n-1)(n-2)(n-3)$ &
$ \scriptstyle \frac{n_c (n_c-1)(n_c-2)(n_c-3)}{n(n-1)(n-2)(n-3)}$ & $\scriptstyle \frac{n_c^2(n_c-1)^2 }{n^2(n-1)^2 }$
& $\scriptstyle \frac{n^2}{n_c^2(n-1)^2}$ & \makecell{$\scriptstyle \hspace{-6mm} \frac{-4 n_t n^3}{n_c^2(n-1)(n_c-1)} -\frac{6n^2}{n_c^2(n-1)} $ \\ $\scriptstyle +\frac{6n^3}{n_c^3(n_c-1)}=O(1) $}
\\ \hline $\scriptstyle i,j \in \textnormal{T}$ $\scriptstyle k,l\in \textnormal{C}$ & $\scriptstyle i=j=k=l$ & $\scriptstyle 2n$ & $\scriptstyle 0$ & $\scriptstyle \frac{n_c n_t}{n^2}$ & $\scriptstyle \frac{n^2 }{n_t n_c}$ & $\scriptstyle -\frac{n^2}{n_t n_c}=O(1)$
\\ or $\scriptstyle i,j \in \textnormal{C}$ $\scriptstyle k,l\in \textnormal{T}$ & \makecell{$\scriptstyle i=j=k,l$ or $\scriptstyle i=j=l,k$ \\ or $\scriptstyle i,j=k=l$ or $\scriptstyle i=k=l,j$} & $\scriptstyle 8n(n-1)$ & $\scriptstyle 0$ & $\scriptstyle \frac{n_c n_t (n_t-1)}{n^2(n-1)}$ & $\scriptstyle-\frac{n^2}{n_t n_c(n-1)}$ & $\scriptstyle -\frac{8n^2}{n_t n_c}=O(1)$
\\ & $\scriptstyle i=j,k=l$ & $\scriptstyle 2 n (n-1)$ & $\frac{n_t n_c}{n^2}$ & $\frac{n_c n_t}{n^2}$ & $\frac{n^2 }{n_t n_c}$ & $\scriptstyle 0$
\\ & $\scriptstyle i=j,k,l$ & $\scriptstyle 4 n(n-1)(n-2)$ & $\frac{n_t n_c(n_t-1)}{n^2(n-1)}$ & $\frac{n_c n_t (n_t-1)}{n^2(n-1)}$ & $ \frac{n^2}{n_t n_c(n-1)}$ & $\scriptstyle 0$
\\ & $\scriptstyle i=k,j=l$ or $\scriptstyle i=l,j=k$ & $\scriptstyle 4n(n-1)$ & $\scriptstyle 0$ & $\frac{n_c n_t (n_t-1) (n_c-1)}{n^2(n-1)^2}$ & $\frac{n^2}{n_t n_c(n-1)^2}$ & $\scriptstyle \frac{4 n^2}{n_t n_c(n-1)}=O(1/n)$
\\ & $\scriptstyle i=k=l,j$ or $\scriptstyle i,j=k=l$ & $\scriptstyle 4n(n-1)$ & $\scriptstyle 0$ & $\frac{n_c n_t (n_c-1)}{n^2(n-1)}$ & $-\frac{n^2}{n_t n_c(n-1)}$ & $\scriptstyle -\frac{4 n^2}{n_t n_c}=O(1)$
\\ & \makecell{$\scriptstyle i=k, j, l$ or $\scriptstyle i=l,j,k$ \\ or $\scriptstyle i,j=k,l$ or $\scriptstyle i,j=l,k$} & $\scriptstyle 8n (n-1)(n-2)$ & $\scriptstyle 0$ & $\frac{n_c n_t (n_c-1) (n_t-1)}{n^2(n-1)^2}$ & $\frac{n^2}{n_t n_c(n-1)^2}$ & $\scriptstyle \frac{8 n^2 (n-2)}{n_t n_c (n-1) }=O(1)$
\\ & $\scriptstyle i,j,k,l$ & $\scriptscriptstyle 2 n(n-1)(n-2)(n-3)$
& $\frac{n_t (n_t-1) n_c (n_c-1)}{n(n-1)(n-2)(n-3)}$ & $\frac{n_c n_t (n_c-1)(n_t-1)}{n^2(n-1)^2}$ & $\frac{n^2}{n_t n_c(n-1)^2}$ & $\scriptstyle \frac{4 n^2 (n-3)}{n_t n_c (n-1) }=O(1)$
\\ \hline $\scriptstyle i,j,k,l \in \textnormal{T}$ & $\scriptstyle i=j=k=l$ & $\scriptstyle n$
& $\frac{n_t }{n}$ & $\frac{n_t ^2}{n^2}$
& $\frac{n^2}{n_t ^2}$ & $\scriptstyle \frac{n^2 n_c}{n_t ^3}=O(1)$
\\ & \makecell{$\scriptstyle i=j=k,l$ or $\scriptstyle i=j=l,k$ \\ or $\scriptstyle i,j=k=l$ or $\scriptstyle i=k=l,j$} & $\scriptstyle 4n(n-1)$ & $\frac{n_t (n_t -1)}{n(n-1)}$ & $\frac{n_t ^2(n_t -1)}{n^2(n-1)}$
& $-\frac{n^2}{n_t ^2(n-1)}$ & $\scriptstyle -\frac{4n^2 n_c}{n_t ^3}=O(1)$
\\ & $\scriptstyle i=j, k=l$ & $\scriptstyle n(n-1)$ & $ \frac{n_t (n_t -1)}{n(n-1)}$ & $\frac{n_t ^2 }{n^2 }$
& $\frac{n^2}{n_t ^2}$ & $\scriptstyle \frac{n_c n^2}{n_t ^3}=O(1)$
\\ & $\scriptstyle i=j,k,l$ or $\scriptstyle i,j,k=l$ & $\scriptstyle 2n(n-1)(n-2)$ & $ \frac{n_t (n_t -1)(n_t -2)}{n(n-1)(n-2)}$ & $\frac{n_t ^2(n_t -1) }{n^2(n-1) }$
& $-\frac{n^2}{n_t ^2(n-1)}$ & $\scriptstyle \frac{4n_c n^2}{ n_t ^3}=O(1)$
\\ & $\scriptstyle i=k,j=l$ or $\scriptstyle i=l, j=k$ & $\scriptstyle 2n(n-1)$ &
$ \frac{n_t (n_t -1)}{n(n-1)}$ & $\frac{n_t ^2(n_t -1)^2 }{n^2(n-1)^2 }$
& $-\frac{n^2}{n_t ^2(n-1)^2}$ & $\scriptstyle -\frac{2 n^2}{n_t ^2}=O(1)$
\\ & \makecell{$\scriptstyle i=k, j, l$ or $\scriptstyle i=l,j,k$ \\ or $\scriptstyle i,j=k,l$ or $\scriptstyle i,j=l,k$} & $\scriptstyle 4n (n-1)(n-2)$
& $ \frac{n_t (n_t -1)(n_t -2)}{n(n-1)(n-2)}$ & $\frac{n_t ^2(n_t -1)^2 }{n^2(n-1)^2 }$
& $\frac{n^2}{n_t ^2(n-1)^2}$ & $\scriptstyle \frac{4 n^2}{n_t ^2}\left(\frac{n(n_t -2)}{n_t (n_t -1)}-\frac{(n-2)}{(n-1)} \right)=O(1)$
\\ & $\scriptstyle i,j,k,l$ & $\scriptscriptstyle n(n-1)(n-2)(n-3)$ &
$ \scriptstyle \frac{n_t (n_t -1)(n_t -2)(n_t -3)}{n(n-1)(n-2)(n-3)}$ & $\scriptstyle \frac{n_t ^2(n_t -1)^2 }{n^2(n-1)^2 }$
& $\scriptstyle \frac{n^2}{n_t ^2(n-1)^2}$ & \makecell{$\scriptstyle \hspace{-6mm} \frac{-4 n_c n^3}{n_t^2(n-1)(n_t-1)} -\frac{6n^2}{n_t^2(n-1)} $ \\ $\scriptstyle +\frac{6n^3}{n_t^3(n_t-1)}=O(1) $}
\\ \hline
\end{tabular}
\caption{Analysis of Condition (ref) for Complete Randomization}
\end{table}
thebibliography\bibitem[\citeauthoryear{Athey and Imbens}{2017}]{atheyimbens} Athey, S., and Imbens, G.W. 2017. The Econometrics of Randomized Experiments. {\it Handbook of Economic Field Experiments}, 1: 73-140.
\bibitem[\citeauthoryear{Arceneaux and Nickerson}{2009}]{arcenick} Arceneaux, Kevin, and David Nickerson. 2009. Modeling uncertainty with clustered data: A comparison of methods, Political Analysis, 17: 177--90.
\bibitem[\citeauthoryear{Aronow and Samii}{2012}]{aronowsamii} Aronow, Peter M. and Cyrus Samii. 2012. Conservative variance estimation for sampling designs with zero pairwise inclusion probabilities. {\it Survey Methodology} 39(1): 231-241.
\bibitem[\citeauthoryear{Aronow and Samii}{2017}]{aronowsamii17} Aronow, Peter M. and Cyrus Samii. 2017. Estimating average causal effects under general interference, with application to a social network experiment. Forthcoming at {\it The Annals of Applied Statistics} .
\bibitem[\citeauthoryear{Aronow and Middleton}{2013}]{aronowmiddleton15} Aronow, Peter M. and Joel A. Middleton. 2013. A class of unbiased estimators of average treatment effect in randomized experiments. {\it Journal of Causal Inference} 1(1): 135-154.
\bibitem[\citeauthoryear{Basse and Feller}{2017}]{bassefeller17} Basse, G. and A. Feller. 2017. Analyzing two-stage experiments in the presence of interference, {\it Journal of the American Statistical Association}, DOI: 10.1080/01621459.2017.1323641
\bibitem[\citeauthoryear{Bloniarz et al.}{2016}]{bloniarzetal} Bloniarz, A., Liu, H., Zhang, C.H., Sekhon, J.S., and Yu, B. 2016. Lasso adjustments of treatment effect estimates in randomized experiments, {\it Proceedings of the National Academy of Sciences}, 113(27): 7383-90.
\bibitem[\citeauthoryear{Campbell and Meyer}{2009}]{campbellmeyer}Campbell, Stephen L., and Carl D. Meyer. 2009. Generalized Inverses of Linear Transformations. {\it https://doi.org/10.1137/1.9780898719048}
\bibitem[\citeauthoryear{Fuller}{2009}]{fuller} Fuller, W.A. 2009. {\it Sampling Statistics.} New Jersey: Wiley.
\bibitem[\citeauthoryear{Fuller and Isaki}{1981}]{fullerisaki} Fuller, W.A. and C.T. Isaki. 1981. {Survey Design Under Superpopulation Models} In: Current Topics in Survey Sampling Eds: Krewski, D. , J.N.K. Rao, R. Platek. New York, Academic Press.
\bibitem[\citeauthoryear{Freedman}{2008a}]{freedman08a} Freedman, D.A. 2008a. On regression adjustments to experimental data. {\it Adv. in Appl. Math.} \textbf{40} 180--193.
\bibitem[\citeauthoryear{Freedman}{2008b}]{freedman08b}Freedman, D.A. 2008b. On regression adjustments in experiments with several treatments. {\it Ann. Appl. Stat.} \textbf{2} 176--196.
\bibitem[\citeauthoryear{Hansen and Bowers}{2009}]{hansenbowers} Hansen, B. and Bowers, J. (2009). Attributing effects to a cluster-randomized get-out-the-vote campaign. {\it J. Amer. Statist. Assoc.} {\bf 104} 873--885.
\bibitem[\citeauthoryear{Holland}{1986}]{holland} Holland, P.W. 1986. Statistics and Causal Inference, \emph{Journal of the American Statistical Association}, vol. 81, no. 396: 945-968.
\bibitem[\citeauthoryear{Horvitz and Thompson}{1952}]{horvitz} Horvitz, D.G. and Thompson, D.J. 1952. A generalization of sampling without replacement from a finite universe. {\it J. Amer. Statist. Assoc.} \textbf{47}: 663-684.
\bibitem[\citeauthoryear{Isaki and Fuller}{1982}]{isakifuller} Isaki, C.T., and W.A. Fuller. 1982. Survey Design Under the Regression Superpopulation Model. \textit{Journal of the American Statistical Association} \textbf{77}(377): 89-96
\bibitem[\citeauthoryear{Li and Ding}{2017}]{liding} Li, Xinran and Ding, Peng. 2017. General Forms of Finite Population Central Limit Theorems with Applications to Causal Inference. {\it Journal of the American Statistical Association} \textbf{112}(520): 1759-1769
\bibitem[\citeauthoryear{Li et al.}{2017}]{lidingrubin} Li, Xinran, Peng Ding, and Donald B. Rubin. 2017. Asymptotic Theory of Rerandomization in Treatment-Control Experiments. \textit{arXiv:1604.00698}
\bibitem[\citeauthoryear{Lin}{2013}]{lin} Lin, Winston. 2013. Agnostic Notes on Regression Adjustments to Experimental Data:
Reexamining Freedman's Critique. {\it Annals of Applied Statistics} \textbf{7}(1): 295-318
\bibitem[\citeauthoryear{Lu}{2016}]{lu} Lu, J. 2016.
Covariate adjustment in randomization-based causal inference for 2k
factorial designs. {\it Statistics & Probability Letters}, \textbf{119}:11–20.
\bibitem[\citeauthoryear{Middleton}{2008}]{middleton08} Middleton, J.A. 2008. Bias of the regression estimator for experiments using clustered random assignment. {\it Stat. Probability Lett.} \textbf{78} 2654--2659.
\bibitem[\citeauthoryear{Middleton}{2018}]{middleton18} Middleton, Joel A. 2018. A unified theory of regression adjustment for design-based inference. {\it arXiv:1803.06011}
\bibitem[\citeauthoryear{Middleton and Aronow}{2015}]{middletonaronow} Middleton, Joel A. and Peter M. Aronow. 2015. Unbiased Estimation of the Average Treatment Effect in Cluster-Randomized Experiments. {\it Statistics, Politics and Policy} \textbf{1}:
\bibitem[\citeauthoryear{Neyman}{1923}]{neyman23} Neyman, Jerzy Splawa, D. M. Dabrowska, and T. P. Speed. [1923.] 1990. On the application of probability theory to agricultural experiments: Essay on principles, section 9. {\it Statistical Science} \textbf{5}: 465--480.
\bibitem[\citeauthoryear{Raj}{1965}]{Raj} Raj, D. 1965. On a method of using multi-auxiliary information in sample surveys. {\it J. Amer. Statist. Assoc.} \textbf{60} 270--277.
\bibitem[\citeauthoryear{Rohde}{1965}]{rhode} Rohde, Charles. 1965. Generalized Inverses of Partitioned Matrices. {\it Journal of the Society for Industrial and Applied Mathematics} \textbf{13}(4): 1033-1035.
\bibitem[\citeauthoryear{Rubin}{1974}]{rubin74} Rubin, Donald. 1974. Estimating causal effects of treatments in randomized and nonrandomized studies. {\it Journal of Educational Psychology} \textbf{66}: 688--701.
\bibitem[\citeauthoryear{Sarndal et al.}{1992}]{sarndaletal} Sarndal, C.-E., B. Swensson, and J. Wretman. 1992. Model Assisted Survey Sampling.
New York: Springer.
\bibitem[\citeauthoryear{Samii and Aronow}{2012}]{samiiaronow} Samii, Cyrus and Peter M. Aronow. 2012. On Equivalencies Between Design-Based and Regression-Based Variance Estimators for Randomized Experiments. {\it Statistics and Probability Letters}. \textbf{82}: 365--370.
\bibitem[\citeauthoryear{Schochet}{2010}]{schochet} Schochet, Peter Z. 2010. Is regression adjustment supported by the Neyman model for causal inference? {\it Journal of Statistical Planning and Inference} \textbf{140}: 246-259.
\bibitem[\citeauthoryear{Sinclair et al.}{2012}]{sinclairetal} Sinclair, B., McConnell, M. and Green, D.P. 2012. Detecting Spillover Effects: Design and Analysis of Multilevel Experiments. {\it American Journal of Political Science} \textbf{56}(4): 1055-1069.
\bibitem[\citeauthoryear{Wood}{2008}]{wood} Wood, John. 2008. On the Covariance Between Related Horvitz-Thompson Estimators. {\it Journal of Official Statistics}. \textbf{24} 53--78.
\bibitem[\citeauthoryear{Zhao et al.}{2017+}]{Zhao} Zhao, A., Ding, P., Mukerjee, R., and Dasgupta, T. 2017+. Randomization-Based Causal Inference From Unbalanced $2^2$ Split-Plot Designs. \textit{Annals of Statistics}, in press.
appendix
Notation Index
tabular[tabular omitted — 2,847 chars of source]
tabular[tabular omitted — 1,489 chars of source]