Design-based Estimation Theory for Complex Experiments

Design-based Estimation Theory for Complex Experiments

Introduction

Randomized experiments have become a standard tool in economic research. Traditionally presented as estimating the average effect of a binary treatment, modern experimental designs have been greatly enriched to capture a variety of economically relevant effects, such as time effects (e.g., athey2022design and roth2021efficient), peer effects (e.g., sacerdote2014experimental), social incentives (e.g., ashraf2018social), and spillover effects (e.g., hudgens2008toward, miguel2004worms and cai2015social). Many such experimental designs involve nonstandard treatment assignment mechanisms and/or interference of treatment status among experimental units according to spatial/network/time proximity.\footnote{By interference, we mean the exposure of one unit to treatment may include other units' assignments. This typically arises when researchers are interested in some spillover effects, e.g. hudgens2008toward and aronow2017estimating. } We refer to these experimental designs are as complex experiments.\footnote{The term complex is borrowed from the survey sampling literature chaudhuri2005survey, where it refers to survey designs that depart from common types of random sampling.} Many researchers analyze experimental data using a regression model with (possibly clustered) robust standard errors. Although such procedures are justifiable for simple experimental designs,\footnote{For example, linear regression models are justifiable in two-arm completely randomized designs (freedman2008regression, lin2013agnostic).} they can be ad hoc when applied to complex experimental designs. It is not clear to what extent the results rely on the modeling assumptions and how to interpret the results when the regression models are thought to be misspecified.

Design-based statistical theory provides a powerful framework for analyzing complex experiments. In the design-based framework, the randomization of treatment assignment is the sole source of statistical randomness. Estimation and inferential theory are formulated on this randomness alone, without reference to any other stochastic model (e.g., sampling from a superpopulation and/or random disturbance terms). This framework has important implications for weighting in the estimation of average treatment effects and for the estimation of standard errors. In simple experiments, the design-based framework provides procedures compatible with current empirical practices for analyzing experimental data with regression models. But in more complex settings, the design-based framework can nevertheless be adapted to provide general-purpose estimation strategies that do not rely on regression models for validity.

Estimation theory in the design-based setting has been investigated for many designs on a case-by-case basis. Many important insights have been derived from studying particular experimental designs, but a design-based estimation theory that can be applied to general experimental designs has not hitherto been developed. A design-based estimation theory with broad applicability is important for practice, as it provides guidance to empirical researchers using novel experimental designs that deviate from well-analyzed cases or simple experimental designs that deviate from the standard ones due to practical limitations and implementation reasons. Such designs appear frequently in economic research.

This paper studies design-based estimation theory for general experimental designs. Our results can be applied to standard designs (e.g., completely randomized designs, clustered randomized designs, and pairwise randomized designs) as well as complex designs where analytical results were not previously available. Under mild regularity assumptions, we provide procedures to consistently and efficiently estimate the average effects of interest and procedures to consistently estimate asymptotic variance bounds.\footnote{In the design-based framework, the asymptotic variance is not generally identified. Starting with splawa1990application, the common solution to the issue of unidentified variances has been to estimate a {\it variance bound}, an identified quantity that is provably greater than the variance. The variance bound formula reduces to the standard (cluster) robust standard errors in simple designs. For example, see lin2013agnostic and schochet2021design.} We also provide a novel scalar measure of experimental complexity which can be linked to the design-based variance of the estimators, enabling researchers to understand the strengths and weaknesses of particular experimental designs. This measure can be used in the designing-stage of the experiment before collecting any outcome data.

\sloppy Building off of recent advances in design-based estimation theory middleton2018unified,middleton2021unifying, the paper makes three main contributions. As the first contribution, we extend the theoretical analysis of many standard estimators to a broader class of experimental designs. Specifically, we analyze a family of design-estimator pairs commonly employed by researchers in practice. We define the class of moment estimators and study their properties with general experimental designs. Special cases of these estimators include the inverse-probability weighted (IPW), Hajek, weighted least squares (WLS), and generalized regression estimators.\footnote{Generalized regression estimators have the same form as doubly-robust estimators in the observational setting, as noted by kang2007demystifying. } We provide conditions for convergence to probability limits and characterize the asymptotic variances for these estimators.\footnote{Refer to Section (ref) for the definition of asymptotics in this setting. } We provide procedures for consistent plug-in variance-bound estimation for general designs under a weak moment assumption.

As a second contribution, we offer new estimators that have desirable asymptotic properties and are applicable with general experimental designs. The new estimators increase estimation precision by having smaller design-based asymptotic variances. The new classes of estimators are based on the class of generalized regression estimators. The first class we consider is the class of standard Quasi-Maximum Likelihood GR estimators (QMLE-GR). This class follows the classical model-assisted estimation strategy in the survey analysis literature sarndal2003model and it is useful when the researcher has a good approximating model for potential outcomes and covariates. However, in terms of asymptotic variances, this strategy is not guaranteed to be superior to the baseline IPW estimator when the model is misspecified. This problem motivates the second class of estimators, the no-harm GR estimators (No-harm-GR). This class of estimators is based on the QMLE estimates but estimates a multiplicative constant in addition. Estimators of this class have an asymptotic variance no worse than that of the baseline IPW estimator. This class of estimators is inspired by the cohen2020no's estimators in a two-arm completely randomized design. The final class is the optimal GR estimators (Opt-GR). This class of estimators leads to the greatest reduction of asymptotic variances when compared with estimators using the same class of parametric models for adjustments. This class of estimators can be traced back to lin2013agnostic, and middleton2018unified studies such estimators for linear models in two-arm experiments. We further consider refinements that combine some of the above approaches, leading to a class of Optimal-Imputed GR estimators (Opt-I GR). We demonstrate the finite sample performances of the proposed estimators using simulated datasets based on an actual network experiment (cai2015social).

As a third contribution, we propose measures of experimental complexity, as a result of the characterization of asymptotic variances. These measures are the largest eigenvalues of the variance-covariance matrices of the inverse probability-weighted treatment assignment indicators. Theoretically, these quantities govern the rate of convergence of moment-type estimators from a design-based point of view. We shall also give a minimax interpretation for such measures: they are the worst-case variance of the IPW estimators when the outcomes are restricted to a unit ball. A collection of such measures provides useful scalar summaries of the relative strengths and weaknesses of an experimental design for measuring different effects of interest. We believe that these measures are useful for researchers to better understand their experimental designs in complex settings and we demonstrate their uses in the simulations.

Literature Review

This paper builds on the insights in middleton2018unified,middleton2021unifying, which proposed the use of matrix spectral theory in the design-based framework. This paper inherits and generalizes the insight. Compared with the previous works, this paper 1) provides a rigorous asymptotic analysis for a large class of estimators (moment-type estimators), 2) considers general asymptotic variance bound estimation under weaker conditions, 3) proposes and analyzes new classes of estimators (QMLE-GR, No-harm-GR and Opt-GR), 4) proposes the measures of complexity and 5) specializes the results to network experiments.

This paper adds to the literature on design-based estimation theory. The survey sampling literature includes a large body of literature on design-based estimation theory (for example, see sarndal2003model and chaudhuri2005survey). Many results in the literature focus on estimating average/total quantities in complex (but not fully general) survey designs and do not consider interference. We consider the case of estimating the contrast of multiple average quantities under general experimental designs and our setup accommodates interference.

We contribute to the literature on estimation theory for the design-based analysis of experiments imbens2015causal. freedman2008b,freedman2008randomization, freedman2008regression, lin2013agnostic, bloniarz2016lasso, wu2018loop, guo2021generalized, cohen2020no and lei2021regression study estimation problems in two-arm completely randomized designs. A collection of papers studies estimation and inferential theory with various experimental designs, for example, middleton2015unbiased, lu2016randomization,li2019rerandomization, roth2021efficient, schochet2021design, negi2021revisiting,athey2022design and gao2023causal. hudgens2008toward, aronow2017estimating, hu2022average and gao2023causal study estimation theory in experiments with interference, and pollmann2020causal studies spatial experiments. aronow2013class considers unbiased difference-type estimation for complex experiments but the paper does not provide any guarantees of variance reduction. Our paper builds on the previous insights and considers the case of general experimental designs, parametric linear and nonlinear models for adjustments, and various strategies for estimating the adjusting models (QMLE, No-harm, Optimal and Optimal-I).

This paper is also related to the literature on variance characterization and variance bound estimation in design-based settings, for example, robins1988confidence,mukerjee2018using, aronow2014sharp,pashley2021insights,de2020level,xu2022. harshaw2021optimized studies the problem of optimizing variance bounds.

Our paper is also related to the literature that analyzes experiment data accounting for variation from both model-based and design-based uncertainties, for example, bugni2018inference, bugni2019inference, bai2021inference, bai2022inference, cytrynbaum2021designing and bugni2022inference. abadie2020sampling provides inferential results for the linear regression model that allows for both design-based and sampling-based uncertainty.

The organization of the paper is as follows. Section (ref) includes model setup, notations and basic assumptions. Section (ref) gives a network experiment example. Section (ref) defines moment-type estimators and provides estimation and variance bound estimation results. Section (ref) includes results on various model-assisted estimation strategies (QMLE, No-harm, Optimal and Optimal-I). Section (ref) provides simulation results based an actual network experiment (cai2015social).

Setup and Notations

We consider a Neyman causal model splawa1990application,imbens2015causal, where one conducts a randomized experiment with $k$ treatment arms on $n$ experimental units. Each unit $i\in\{1,...,n\}$ is associated with a $k$-vector of nonrandom potential outcomes:

Each unit $i$ is randomly assigned to one of the $k$ treatment arms. We denote the random vector of assignment indicators by

where ${D}_{ai}=1$ means that the unit $i$ is assigned to the treatment arm $a$ and ${D}_{ai}=0$ otherwise. For unit $i$, the observed outcome is generated according to

One may also observe for each unit $i$ an additional $p$-dimensional row vector of pretreatment covariates $x_i=(x_{1i},x_{2i},...,x_{pi})\in\mathbb{R}^{p}$. In this paper, we assume that the dimension of the covariates does not change with the sample size. We stack the covariate vectors vertically to create a matrix $\mathbf{x}\in\mathbb{R}^{n\times p}$. The observed data for unit $i$ can then be represented as $(Y_i^{{\scriptscriptstyle{\textnormal{obs}}}},{D}_{1i}, {D}_{2i},...,{D}_{ki},x_i)\in\mathbb{R}^{1+k+p}$. Let $\pi_{ai}=\text{\textnormal{E}}[D_{ai}]$ denote the probability of assignment of unit $i$ to treatment arm $a$. We shall hereafter assume that $\pi_{ai}$ is positive for all treatment arm $a$ and unit $i$, unless stated otherwise.

The parameters of interest are contrasts (linear combination) between the group-specific means of the potential outcomes. For example, in a two-arm experiment, the parameter of interest could be the average treatment effect (ATE) between the treated group and control group and it is defined as $\frac{1}{n} \sum_i \left(y_{i}(2)-y_{i}(1)\right)$.

Notation

Let $y^1$, $y^2$, ..., $y^k$ represent column $n$-vectors of potential outcomes associated with each of the arms, with the $i$th element of each vector corresponding to the $i$th unit. Thus, $y^a=[y_{i}(a)]_{i=1}^n=\left(y_{1}(a),y_{2}(a),...,y_{n}(a)\right)'\in\mathbb{R}^n$. We stack these vectors vertically to create a column vector $y$ with length $kn$:

Let $1_{\scriptscriptstyle n}$ be a column $n$-vector of ones. A $kn \times k$ intercept matrix is defined as {

} Entries left blank are equal to 0.\footnote{Formally, this matrix is defined as $\mathbf{1}=[a_{st}]_{s=1,..,kn}^{t=1,...,k}$, where $a_{st}=1$ if $(t-1)\leq \frac{s}{n}\leq t$ and 0 otherwise. } A k-vector of the average potential outcomes of the arms can then be written as $\mu_n=\frac{1}{n} \mathbf{1}' y$. From here on, we denote the $k$-vector average potential outcomes as $\mu_n$ and estimators as $\hat{\nu}_n^{\textnormal{(type)}}\in\mathbb{R}^k$. The superscript denotes the type of estimator. For example, an inverse-probability weighted estimator for the average potential outcomes will be denoted as $\hat{\nu}_n^{{\scriptscriptstyle{\textrm{IPW}}}}$. \\ Next, define an $n \times n$ diagonal matrix with $n$ assignment indicators for treatment arm 1 on the diagonal and 0 otherwise,

and define ${D}^2$, ${D}^3$, $\hdots$, ${D}^k$ analogously. Arrange these matrices to create a diagonal $kn \times kn$ matrix {

} The $kn\times kn$ diagonal matrix of assignment probabilities is written as $\boldsymbol{\pi}=\text{\textnormal{E}}[{D}]$. \\ For the purpose of covariate adjustments, we also define the $kn \times (k+p) $ matrix, {

} which augments the intercept matrix $\mathbf{1}$ with covariates.

Let $c\in\mathbb{R}^k$ denote an arbitrary column contrast vector such that the parameter of interest can be written as $\frac{1}{n}c'\mathbf{1}' y$. For example, in a two-arm experiment, the ATE is defined by choosing $c=\left(-1 ,1\right)'$ and $\frac{1}{n}c'\mathbf{1}' y=\frac{1}{n} \sum_i \left(y_{i}(2)-y_{i}(1)\right)$. We write the parameter of interest associated with a contrast vector $c$ as $\mu_n^c= \frac{1}{n}c'\mathbf{1}' y$.

To conclude, in this notation, we say researchers observe the assignment ${D}$, outcomes ${D} y$, and a matrix of $p$ pretreatment covariates $\mathbf{x}\in\mathbb{R}^{n\times p}$. In a randomized experiment, $\boldsymbol{\pi}$ is known, or can be approximated to arbitrary precision by repeating the randomization procedure and collecting draws fattorini.

We let $I_{k}$ denote the identity matrix of dimension $k\times k$ and $\mathbf{0}_{k}$ a zero matrix of dimension $k\times k$. $ 1_{\scriptscriptstyle {k}}$ denotes a column $k$-vector of 1s. We define $[k]=\{1,...,k\}$ and $[k_1,k_2]=\{k_1,...,k_2\}$, for arbitrary positive integers $k$, $k_1$ and $k_2$. A list of mathematical objects, operators, and quantities used in this paper is included in Appendix (ref).

Asymptotic Schemes

All results in this paper are asymptotic. We consider a nested sequence of increasing populations, $\{U_n\}_{n}$, where the index $n$ indicates the size of the population under study. Each unit in the population is characterized by its fixed potential outcomes and pretreatment covariates. The potential outcomes and the pretreatment covariates are fixed and the population grows deterministically. $\{U_n\}_{n\geq 1}$ are nested: $U_1\subset U_2\subset ...\subset U_n...$. Each finite population $U_n$ has an associated experimental design and a realized randomization. Although the populations form a nested sequence, the sequences of realized assignments do not. This asymptotic scheme is widely used in the literature isaki1982survey,aronow2014sharp,li2017general.

We work in the finite-population (design-based) framework imbens2015causal, where the potential outcomes $y$ and pretreatment covariates $\mathbf{x}$ are considered fixed parameters, and the only source of randomness in our model is from the random treatment assignments $\{D_{ai}\}_{a\in[k],i\in[n]}$. \\ In general, the true parameter values are quantities that change with the sample size $n$. We will write (finite) population quantities with a subindex $n$. For example, the average potential outcomes will be denoted as $\mu_n=\frac{1}{n}\mathbf{1}'y$. We call $\hat{\nu}_n$ a consistent estimator for $\nu_n$ if $\hat{\nu}_n-\nu_n=o_p(1)$, where the stochasticity is generated by the experimental design. With an abuse of language, we refer to $\nu_n$ as the probability limit of $\hat{\nu}_n$.

We state two assumptions for data moments, which are needed for the convergence of estimators. Recall that $k$ denotes the number of treatment arms and $p$ denotes the number of pretreatment covariates.

For our analysis of weighted least squares (WLS) estimators below, we require the design matrix to be invertible for large $n$.

An example: network experiments

This section provides a concrete example to illustrate the setup described above. We consider the network experiments proposed in aronow2017estimating. Components of the experimental design include:

For one experiment, researchers randomly draw an assignment vector $Z$ and observe the scalar outcomes $\{Y_i(Z)\}_{i=1}^n$. Notice that up to this stage the outcome for unit $i$ has depended on the entire assignment vector. Consistent estimation under unrestricted interference is deemed virtually impossible savje2021average. One strategy to alleviate the problem is to restrict the interference patterns using exposure mappings, under the assumption that the potential outcomes are correctly specified according to the exposure value:

Assumption (ref) implies that the exposures are "effective treatments" as defined in manski2013identification. The assumption that $\Delta_n$ equals a finite set $\Delta$ is a typical assumption made in the literature. We note that this experimental setup is very general, and it can be generalized to other settings in which the exposure mappings are not necessarily mediated by a network.\\ Enumerate the element in $\Delta$ as $\{1,...,k\}$. With Assumption (ref), one can write the potential outcomes associated with unit $i$ as $\left(y_{i}(1),...,y_{i}(k)\right)\in\mathbb{R}^k$. The assignment vector associated with unit $i$ can be written as $\left({D}_{1i},...,{D}_{ki}\right)$. This maps the problem back to our framework. To make the assumptions concrete, we provide an example from Section 9 of aronow2017estimating.

The theories examined in this paper can be viewed as offering strategies for analyzing and planning such (though not limited to) experiments.

Estimation and Variance Estimation

In this section, we study the problem of estimating average potential outcomes, $\mu_n=\frac{1}{n}\mathbf{1}'y\in\mathbb{R}^k$.

Section (ref) introduces the class of moment-type estimators that allows the simultaneous analysis for many commonly-used estimators. It nests the class of linear estimators introduced in mukerjee2018using and middleton2021unifying, which includes the IPW, Hajek (HJ), Weighted Least Square (WLS), Completely Imputed (CI), Missing Imputed (MI) and Generalized Regression (GR) estimators. We establish the convergence rate of the moment-type estimators, and give conditions under which the WLS is a consistent estimator of the average potential outcomes.

In Section (ref), we study asymptotic variance characterization, bounding, and variance bound estimation for the moment-type estimators. We provide a simple and general asymptotic variance formula. We highlight that the asymptotic variances can be written in the bilinear form $\frac{1}{n}z'\mathbf{\Omega} z$, where $z$ reflects the choice of estimators and the matrix $\mathbf{\Omega}$, introduced in Definition (ref), reflects the experimental design. This forms the basis for variance bounding and variance bound estimation. We discuss variance bounding and provide a plug-in variance bound estimator. We briefly discuss inference in Section (ref).

Appendix (ref) formulates the lower-level conditions for the network experiments presented in Section (ref).

In Section (ref), we propose using the largest singular value $\sigma_{\max}\left(\mathbf{\Omega}\right)$of the matrix $\mathbf{\Omega}$ as an input for experimental designs. The value can be interpreted as the worst-case variance of the IPW estimator when the outcomes are subject to a moment condition.

We remind readers that $a$ is an index for treatment arms, $i$ is an index for experiment units, $k$ is the number of treatment arms, and $n$ is the number of experiment units.

Moment-type estimator

We first define moment-type estimators. We will give several examples of moment-type estimators after the definition.

We associate a probability target for each moment-type estimator.

In short, the probability limit of a moment-type estimator is obtained by replacing estimated moments with their finite population counterparts. For notational simplicity, we hereafter write $\widehat{m}_n=\allowdisplaybreaks(\widehat{m}_n^1,\widehat{m}_n^2,...,m_n^{S_1+S_2})$, $m_n=\allowdisplaybreaks(m_n^1,m_n^2,...,m_n^{S_1+S_2})$. We shall write $\widetilde{m}_n=\allowdisplaybreaks (\widetilde{m}_n^1,\widetilde{m}_n^2,...,\widetilde{m}_n^{S_1},m_n^{S_1+1},...,m_n^{S_1+S_2})$ for some arbitrary $\widetilde{m}_n^1, ...,\widetilde{m}_n^{S_1}$.

The moment-type estimators nest many commonly-used estimators as special cases.

All estimators listed above are moment-type estimators. For example, the IPW estimator can be written as:

which are a mapping of unknown moments of potential outcomes.\\ The WLS estimator can be represented as

where $\{\widehat{\beta}_{a}\}_{a\in[k]}$ are the realized estimators of the arm-specific coefficients from the weighted regression:

The WLS estimator is a mapping of terms such as $\frac{1}{n}\sum_{i=1}^n {D}_{ai}x_{ki}y_{ai}\mathbf{\omega}_{ai}$ or $\frac{1}{n}\sum_{i=1}^n {D}_{ai}x_{ki}x_{si}\mathbf{\omega}_{ai}$ for some $k,s\in[p]$.

The GR estimator takes the doubly-robust influence function form and is commonly considered in the literature. Besides the estimated moments, the GR estimator also involves known moments of the form $\frac{1}{n}\sum_{i=1}^n x_{ki}$, $k\in[p]$.

Besides the given examples, the moment-type estimator also includes Hajek, completely-imputed (CI), and missing-imputed (MI) estimators. These estimators are routinely considered in the survey sampling literature, and we include their definitions in Appendix (ref).\footnote{Missing imputed estimators and completely-imputed estimators are considered in isaki1982survey. CI, MI, and GR estimators have a long history in the survey sampling literature (see, for example, brewer1979class, wright1983finite, sarndal1984cosmetic, and chaudhuri2005survey). Missing imputed estimators have been recently studied by guo2021generalized in nonlinear adjustment problems. }

We now show convergence of the estimator $\widehat{\nu}_n$ to its probability target $\nu_n$. We note that $\nu_n$ is not necessarily the average potential outcomes, and we will specify conditions under which the WLS is a consistent estimator for the average potential outcomes.\footnote{IPW, GR and Hajek estimators have average potential outcomes as their probability target and hence are consistent estimators. Conditions under which the CI and MI estimators are consistent estimators of the average potential outcomes are given in Section (ref). } To characterize the rate of convergence, we introduce the first-order design matrix. This matrix is first introduced in middleton2021unifying and it is an important conceptual object encoding information about the experimental design.

We give an example of the first-order design matrix with two treatment arms $\left(a=1,2\right)$ and two units $\left(i=1,2\right)$: {

} where the subscripts follow the convention ${D}_{ai}$, where $a$ indexes treatment arms and $i$ indexes units, and $\pi_{ai} = \text{\textnormal{E}}[{D}_{ai}]$. It should be clear that the first-order design matrix depends on both first-order and second-order assignment probabilities.

We use $\sigma_{\max}(A)$ to denote the largest singular value of a matrix $A$. We use $\|v\|_2$ to denote the Euclidean ($l_2$) norm when $v$ is a vector and the Frobenius norm when $v$ is a matrix. We show that $\sigma_{\max}(A)$ upper-bounds the statistical convergence rate for the moment-type estimators under weak regularity conditions.

Theorem (ref) implies the following corollary for the IPW, WLS, and GR estimators.

We conclude this section with a lemma establishing the consistency of the WLS estimators for the average potential outcomes. This result implies that the inverse-probability weighted WLS estimator with centered covariates is consistent for estimating the average potential outcomes in general. Results for CI and MI estimators are included in Section (ref).

Asymptotic Variances: Characterization, bounding and estimation

To characterize the asymptotic variance, we strengthen Assumption (ref) to use a linearization argument. In addition to using the notation $m_n$ and $\widehat{m}_n$ from Definition (ref), we further define the column vectors $\widetilde{m}_{n,r}=(\widetilde{m}_n^1,...,\widetilde{m}_n^{l_1})'\in\mathbb{R}^{S_1}$ and $m_{n,r}=(m_n^1,...,m_n^{l_1})'\in\mathbb{R}^{S_1}$, as well as $\widetilde{m}_{n}=(\widetilde{m}_n^1,...,\widetilde{m}_n^{l_1},m_n^{S_1+1},...,m_n^{S_1+s_2})'\in\mathbb{R}^{S_1+S_2}$\footnote{The subscript r stands for "random moments".}

For a moment-type estimator $\widehat{\nu}_n=F(\widehat{m}_n)$, we denote the linearized version of $\widehat{\nu}_n$ as $\widehat{\nu}_n^{\scriptscriptstyle{\textnormal{L}}}=F(m_n)+dF_{m_n}(\widehat{m}_{n,r}-m_{n,r})$, where $\widehat{m}_{n,r}=(\widehat{m}_n^1,...,\widehat{m}_n^{S_1})'\in\mathbb{R}^{S_1}$. Let $\operatorname{diag}()$ be the operator that maps a length-n vector to an n-by-n diagonal matrix.

This theorem is the key result of the section. It is shown that many linearized estimators have a variance that can be written as a "quadratic" form.\footnote{A typical definition of quadratic forms is the expression $v'\mathbf{A}v$, where $v$ is a column vector and $\mathbf{A}$ a square matrix. In the theorem below, our $v$ is a matrix instead of a column vector, so we are in a sense mis-using the name. If we are interested in a scalar parameter specified by a contrast vector $c$, then the variance $\text{\textnormal{Var}}(\sum_{a=1}^k c_a\widehat{\mu}_a^{{\scriptscriptstyle{\textnormal{L}}}})$ is a quadratic form in the standard sense. } Notice $\mathbf{\Omega}$ depends only on the information of experimental designs, which is available to researchers. The matrix $z$ can depend on potential outcomes, the experimental design, and/or covariates, and it needs to be recovered from the data.

middleton2021unifying formally shows that the IPW estimator has an asymptotic variance that can be written as a bilinear form. Our results generalize middleton2021unifying's intuition regarding asymptotic variances with a rigorous proof. With an abuse of language, we refer to $n\text{\textnormal{Var}}(\widehat{\nu}_n^{{\scriptscriptstyle{\textnormal{L}}}})=\frac{1}{n}z'\mathbf{\Omega} z$ as the asymptotic variance of the estimator $\widehat{\nu}_n$. We specialize Theorem (ref) for estimators introduced in Section (ref).

Now we turn to the subject of variance (bound) estimation. In the design-based framework, the true asymptotic variance is not identified nor consistently estimable. One can read off the lack-of-identification problem from entries in the first-order design matrix $\mathbf{\Omega}$: some entries have the value $-1$, and this happens if $y_{ai}$ and $y_{bj}$ can never be simultaneously observed across all realized assignments and we have $\text{\textnormal{E}}[\left({D}_{ai}{D}_{bj}\right)/\left(\pi_{ai}\pi_{bj}\right)]-\text{\textnormal{E}}[{D}_{ai}/\pi_{ai}]\text{\textnormal{E}}[{D}_{bj}\pi_{bj}]=0-1=-1$. To construct a variance bound, one needs to find a variance bound matrix $\widetilde{\mathbf{\Omega}}$ that dominates $\mathbf{\Omega}$ in the positive semidefinite sense. In order for $\widetilde{\mathbf{\Omega}}$ to be identified, $\widetilde{\mathbf{\Omega}}$ must take the value $0$ at entries that are $-1$ in $\mathbf{\Omega}$. The concepts are formalized in mukerjee2018using and middleton2021unifying.

Entries in $\text{\textnormal{I}}\left(\mathbf{\Omega}=-1\right)$ are indications that the associated terms in the variance quadratics are impossible to observe. $\text{\textnormal{I}} (\mathbf{\Omega}=-1) \leq \text{\textnormal{I}}(\widetilde{\mathbf{\Omega}}= 0)$ are indications that the associated terms are not used in the variance bound estimation and thus the variance bound (of an IPW estimator) is identified and can be estimated. We note that it is not always necessary to bound the entire matrix $\mathbf{\Omega}$. Depending on the parameter of interest, one may only need to bound a principal submatrix of $\mathbf{\Omega}$. This happens, for example, if the researcher is only interested in comparing two arms of a multi-arm experiment. From now on, for an estimator $\widehat{\nu}_n$ with an asymptotic variance $\frac{1}{n}z'\mathbf{\Omega} z$, we refer to the estimator's asymptotic variance bound as $\frac{1}{n}z'\widetilde{\mathbf{\Omega}} z$.\footnote{ To demonstrate how this definition maps to a typical case, consider a two-arm completely randomized experiment with $n_t$ units in the treatment group and $n_c$ units in the control group, with $n=n_t+n_c$. Define the rescaled demeaning matrix $\mathbf{A}_n=\frac{n}{n-1}\left(I_{n}-\frac{1}{n} 1_{\scriptscriptstyle {n}} 1_{\scriptscriptstyle {n}}'\right)\in\mathbb{R}^{n\times n}$. The first-order design matrix for the design is

The standard Neyman bound (e.g. see imbens2015causal) matrix can be written as

Note that no entries in $\widetilde{\mathbf{\Omega}}^{{\scriptscriptstyle{\textnormal{N}}}}$ take on the value $-1$ and the added matrix is a positive semidefinite matrix. Thus, the Neyman bound matrix is an identified variance-bound matrix by the definition.}

For our purpose, we need a variance-bound matrix that can suit general experimental designs. The Aronow-Samii variance bound (aronow2017estimating) is a general variance bound and applicable to arbitrary designs. We refer readers to harshaw2021optimized for a general variance-bounding technique.

With an identified variance-bound matrix $\widetilde{\mathbf{\Omega}}$ (not necessarily the AS bound), we turn to variance bound estimation. We provide regularity conditions for consistent plug-in variance bound estimation. We first define the second-order design tensor.\footnote{The tensor is first introduced in middleton2021unifying.} The operator norm of this object, which we define below, determines the rate of convergence for the variance-bound estimator. We use $\otimes$ to denote the tensor product of two matrices, which results in an fourth-order tensor. For any two matrices $A=[a_{ij}]\in\mathbb{R}^{n_1\times n_2}$ and $B=[b_{ij}]\in\mathbb{R}^{n_3\times n_4}$, we denote

Next, define an inverse probability weighted version of the variance bound matrix, $\widetilde{ \mathbf{\Omega} }$, as

where $\mathbf{p}$ is defined in Definition (ref) and $/$ denotes elementwise division with division by zero defined to be zero. With the matrix, an unbiased estimator of the variance bound $\widetilde{\text{\textnormal{Var}}}\left(\widehat{\nu}_n^{{\scriptscriptstyle{\textnormal{L}}}}\right) = \frac{1}{n^2}z'\widetilde{\mathbf{\Omega}}z$ can be written as

had $z$ is known. Since $z$ generally involves unknown quantities, an appeal to the plug-in principle suggests the use of

where $\widehat{z}$ has the same form as $z$ but with unknown quantities replaced by their estimators.\footnote{ middleton2021unifying suggests the use of plug-in estimators but does not offer a formal justification.}\\ Specializing to the IPW, WLS and GR estimators, the variance bound estimators have the form:

where ${D}\widehat{\epsilon}={D}\left(y-X\widehat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}\right)$ and we note that ${D}\operatorname{diag}\left(\widehat{\epsilon}\right)$ and ${D}\operatorname{diag}\left(y\right)$ do not involve unobserved quantities.

The following theorem considers the problem of consistent plug-in variance bound estimation. We introduce a few tensor notations used below.\footnote{We only introduce the minimally necessary notation here. For complete notation, refer to Appendix (ref).} A real fourth order tensor $\mathbf{A}= \left ({a_{ \scriptscriptstyle i_{\scriptscriptstyle 1}...i_{\scriptscriptstyle 4}}}\right)\in\mathbb{R}^{n_1\times ...\times n_4} $ is a multi-array of entries, where $i_j=1,...,n_j$ for $j=1,...,4$. When $n=n_1=...=n_4$, $\mathbf{A}$ is called a fourth-order $n$-dimensional tensor. For a fourth-order n-dimensional tensor $\mathbf{A}$, we use the symbol $\sigma_{\max}(\mathbf{A})$ to denote the optimal value of the following optimization problem:

This quantity is defined in lim2005singular as a generalization of matrix singular values to tensors. Note that the arguments are constrained to be in the $l_4$ ball instead of the $l_2$ ball in $\mathbb{R}^n$.

Recall that $\mathbf{S}$ is defined in Definition (ref). Let $\|\cdot\|_2$ denote the Frobenius norm if applied to a matrix and the $l_2$ vector norm if applied to a vector. Define $\circ$ to be the entrywise multiplication of two tensors. The following theorem provides the convergence rate for the plug-in variance bound estimator.

We specialize the theorem for the IPW, WLS and GR estimators.

Inference

For hypothesis testing and confidence interval construction, standard theoretical and empirical practice is to establish a central limit theorem for the linearized estimators. For example, let $c\in \mathbb{R}^k$ be a contrast vector, one usually wants to prove $\left(c'\widehat{\nu}_n^{{\scriptscriptstyle{\textnormal{L}}}}-c'\nu_n\right)/\sqrt{c'\text{\textnormal{Var}}\left(\widehat{\nu}_n^{{\scriptscriptstyle{\textnormal{L}}}}\right)c}\overset{d}{\to} N(0,1)$, where $N(0,1)$ is a standard random variable with mean $0$ and variance 1. This justifies the usual inferential procedure for using a normal quantile for testing and confidence interval construction.

Unfortunately, we are not aware of a CLT general enough to be applied with the assumptions introduced so far. It is difficult to prove a CLT in our setting because the assignment random variables are allowed to be almost arbitrarily correlated. In fact, Proposition 11 in savje2021average has shown that the Chevbyshev's inequality is asymptotically sharp when there is strong interference.

For the network experiments discussed in Section (ref), the CLT can be justified when i) there are many disjoint components (e.g. villages and classrooms) in a large network and assignments are independent across different components but can be arbitrarily correlated within; ii) the network can be fully connected but the network has a bounded degree and assignments are independent among units.

When a CLT is not applicable or cannot be rigorously justified, one can still use other tail bounds, for example, the Chebyshev inequality, for inference, although this may entail a loss of efficiency when a CLT indeed approximately holds. The following theorem provides a justification for the plug-in inference. It says that if a critical value can be used for infeasible inference with the linearized estimator, the plug-in inference with the moment-type estimators and variance-bound estimator will remain asymptotically valid under regularity conditions.

In particular, the theorem holds for $\widehat{\nu}_n^{{\scriptscriptstyle{\textrm{IPW}}}}$, $\widehat{\nu}_n^{{\scriptscriptstyle{\textnormal{WLS}}}}$ and $\widehat{\nu}_n^{{\scriptscriptstyle{\textnormal{GR}}}}$ under the premises of Corollary (ref) and if there exist a positive integer $n_0$ and a positive constant $\epsilon$ such that $n\text{\textnormal{Var}}(\nu^{{\scriptscriptstyle{\textnormal{L}}}}_n)\succeq \epsilonI_k$ for all $n\geq n_0$.

$\sigma_{\max}\left(\mathbf{\Omega}\right)$ as an input for designing experiments

The scalar value $\sigma_{\max}\left(\mathbf{\Omega}\right)$ is the largest singular value of the first-order design matrix $\mathbf{\Omega}$ and it appears in the convergence rate calculation in ((ref)). A smaller $\sigma_{\max}\left(\mathbf{\Omega}\right)$ generically implies a better estimation quality for moment-type estimators. We now show that $\sigma_{\max}\left(\mathbf{\Omega}\right)$ can be interpreted as the worst-case variance for the IPW estimator. Recall the variance representation of the IPW estimator:

With a (column) contract vector $c\in\mathbb{R}^k$, the variance of the IPW estimator $c'\widehat{\mu}_n^{{\scriptscriptstyle{\textrm{IPW}}}}$ is:

Since $\mathbf{\Omega}$ is a semidefinite variance-covariance matrix, the largest singular value is the largest eigenvalue. The variational definition of the eigenvalues implies that:

Hence $\sigma_{\max}\left(\mathbf{\Omega}\right)$ can be interpreted as the worst-case variance of the IPW estimator when the outcome is restricted such that $\frac{1}{n}\sum_{i=1}^n\sum_{a=1}^k \left(c_{a}y_{ai}\right)^2\leq 1$.

The quantity introduced above may be too pessimistic because it measures the worst-case scenario with respect to all possible varying contrast vectors $c$ and the potential outcome vector $y$. For a fixed contrast vector $c$, we can define a contrast-specific first-order design matrix:

where $\mathbf{\Omega}_{aa'ii'}=\text{\textnormal{Cov}}\left({D}_{ai}/\pi_{ai},{D}_{a'i'}/\pi_{a'i'}\right)$ is the entry in $\mathbf{\Omega}$ associated with unit $i$ and treatment arm $a$, and unit $i'$ and arm $a'$. It can be shown that:

as a result of the algebraic identity:

and again apply the variational definition of eigenvalues as in ((ref)).

Hence, when we make no additional assumptions other than a scale restriction on the potential outcome vector $y$, the worst-case variance of the IPW estimator with a specific contrast vector $c$ is proportional to $\sigma_{\max}\left(\mathbf{\Omega}^c\right)$. In this sense, if researchers perceive themselves as lacking a reliable prior on the potential outcomes, an experimental design with a smaller $\sigma_{\max}\left(\mathbf{\Omega}^c\right)$ is preferred to one with a larger $\sigma_{\max}\left(\mathbf{\Omega}^c\right)$, when the goal is to measure the average effects defined by the the contrast vector $c$.\footnote{If the researchers have a prior on the potential outcomes, they can instead consider the anticipated asymptotic variance (isaki1982survey), that is, $\frac{1}{n}\mathcal{E}[ z'\mathbf{\Omega} z]$, where $\mathcal{E}$ is the expectation over the randomness with respect to the prior.}

The following numerical example from middleton2021unifying is instructive to conceptualize the measure $\sigma_{\max}\left(\mathbf{\Omega}^c\right)$. We compare $\sigma_{\max}\left(\mathbf{\Omega}^c\right)$ for two treatment-control experiment designs with four units. We compare a completely randomized design where two units are assigned to the treatment group and two to the control group, with a pairwise randomized design where the first two units form a pair and the last two units form a pair. The goal is to measure the average treatment effect defined by the contrast vector $c=(-1,1)$.

The $\mathbf{\Omega}^c$'s for the two designs are: {

}

We have the value 2.667 for $\sigma_{\max}\left(\mathbf{\Omega}^{c,\textrm{cr}}\right)$ and 4 for $\sigma_{\max}\left(\mathbf{\Omega}^{c,\textrm{pair}}\right)$. Hence in terms of robustness (worst-case variance), the completely randomized design is preferable to the pairwise randomized design. This comparison reflects the well-known fact that pairwise randomization is less efficient when potential outcomes are negatively correlated within pairs, as it depends on the researcher's prior belief that paired units have similar outcomes.

We note that this comparison should not be interpreted as a negative result for the matched-pair design. Instead, our primary goal is to underscore the trade-off between robustness and efficiency in experimental designs. In complex experimental settings where prior information on potential outcomes is limited or unavailable, we hope that the proposed metrics $\sigma_{\max}\left(\mathbf{\Omega}\right)$ and $\sigma_{\max}\left(\mathbf{\Omega}^c\right)$ are informative for guiding design choices. We demonstrate the use of these measures in our simulation section in Section (ref).\footnote{There is a long literature on the minimax design of experiments wu1981robustness,li1983minimaxity,kallus2018optimal, harshaw2019balancing, bai2023randomize,basse2023minimax,ni2023design. Our main contribution is to introduce a worst-case variance measure applicable to general experimental designs. Developing a general minimax experimental design is beyond the scope of the present paper.}

Model-Assisted Estimators and Optimality

It can be shown that many commonly used estimators, including properly specified WLS estimators, are GR estimators.\footnote{We include a proof in Appendix (ref).} GR estimators have the doubly-robust form

where $f(X,\widehat{\theta}_n)=X'\widehat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_n$. One intuition reflected in this estimator is that if $X'\widehat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_n$ "predicts" the potential outcomes well, GR estimators may be expected to be more precise relative to the baseline IPW estimator in terms of (asymptotic) variances. This viewpoint motivates the model-assisted estimation strategy in the survey sampling literature (sarndal2003model). However, if a model-assisted estimator is not constructed carefully, it may reduce asymptotic precision when compared with the IPW estimator freedman2008regression,lin2013agnostic.

In this section, we study the problem of model-assisted estimation strategies in general experimental designs using GR estimators. We focus on parametric models in this paper and leave the study of nonparametric or high-dimensional models for future work. We examine three classes of estimators and an additional hybrid class of estimators.\footnote{Each class has some precedence in the literature for some particular experimental designs, and we extend them to general experimental designs and consider a larger class of models.}

The first class we consider is the standard Quasi-Maximum Likelihood GR estimators (QMLE-GR). This class includes the most commonly used models in the literature, such as linear, probit, and logit models. However, in terms of asymptotic variances, this estimation strategy is not guaranteed to be superior to a vanilla IPW estimator freedman2008b,lin2013agnostic. This problem motivates the second class of estimators, the no-harm GR estimators (No-harm-GR). This class is based on the QMLE estimates but estimates a multiplicative constant in addition. This class of estimators provably yields an asymptotic variance no worse than that of the baseline IPW estimator. The final class is the optimal GR estimators (Opt-GR). This class of estimators is optimal in that it achieves the greatest reduction in asymptotic variance among all GR estimators that use the same class of imputation functions (but differ in their parameter values).

Although Opt-GR estimators have strong asymptotic theoretical guarantees, we identify two limitations. One is that they tend to be unstable when $\sigma_{\max}\left(\mathbf{\Omega}\right)$ is large, resulting in poor finite-sample performances. Another drawback is that they require solving nonconvex optimization problems, which may limit their compatibility with, for example, modern machine learning toolboxes. As a result, we propose a hybrid class of estimators called optimal-imputed GR estimators (Opt-I GR). This estimator combines the insights from the No-harm GR and Opt-GR estimators: they are Opt-GR estimators with linear imputation functoins, but instead of using the full set of covariates, the linear models use a single covariate that is imputed by a QMLE model. We find that this class of estimators has good finite-sample performance in our simulations.

We remind readers that $k$ denotes the number of treatment arms and $n$ denotes the number of experiment units. Treatment arms are generically indexed by $a\in[k]$, and units are generically indexed by $i\in[n]$. In many cases, we write $\sum_{a=1}^k\sum_{i=1}^n$ as $\sum_{a,i}$ for simplicity.

Our GR estimators are constructed using the imputation functions $f^a(\cdot,\theta): \mathcal{X}\to \mathbb{R},a\in[k]$ indexed by a finite-dimensional coefficient $\theta$, and $\mathcal{X}$ is the domain of the covariates. Different treatment arms can have their own imputation functions. All our GR estimators are constructed as follows:

Different choices of the criterion $\mathcal{L}(\theta)$ give rise to different estimators. This section examines these choices and their implications.

We define additional notation for later use. Define the column vector $f^a(\theta)=\left(f^a(x_1,\theta),...,f^a(x_n,\theta)\right)\allowbreak\in\mathbb{R}^n$ and column vector $f(\theta)=\left(f^1(\theta)',f^2(\theta)',...,f^k(\theta)'\right)'\in\mathbb{R}^{kn}$. In other words, $f^a(\theta)$ includes the imputed potential outcomes for all units in arm $a$ with parameter $\theta$, and $f(\theta)$ include the imputed outcomes of all arms. We first examine the QMLE-GR estimators and then use the results to motivate the no-harm GR and optimal-coefficient GR estimators.

QMLE-GR estimators

For the QMLE-GR estimators, the population and sample criteria are defined as follows:

The functions $g^a(\cdot,\theta)$ are measures of loss (e.g., squared loss or log likelihood), $\omega_{ai}'s$ are nonnegative weights. For a parameter of interest $\mu_{n,c}=c'\frac{1}{n}\mathbf{1}'y$, we construct the QMLE-GR estimator $\widehat{\mu}_{n,c}^{{\scriptscriptstyle{\textnormal{QMLE}}}}$ as in Algorithm (ref).

\newline The following theorem collects the standard asymptotic results for the QMLE-GR estimators. It shows that, under regularity assumptions, i) the estimator $\widehat{\mu}_{n,c}^{{\scriptscriptstyle{\textnormal{QMLE}}}}$ is equivalent to its linearized counterpart, ii) the plug-in asymptotic variance bound estimator is consistent, iii) the estimator and variance bound estimator can be used for inference. The regularity assumptions (Assumptions (ref) and (ref)) are standard newey1994large, and we include them in Appendix (ref). For simplicity of notations, we shall focus on the $\sqrt{n}$-case for all our estimators below.

Define $\theta_n\equiv\arg\min_{\theta\in\Theta} \mathcal{L}_n(\theta)$, where $\mathcal{L}_n(\theta)$ is defined in ((ref)). Further, define the infeasible estimator:

for $a\in[ k]$ and $\widehat{\mu}_{n,c}^{{\scriptscriptstyle{\textnormal{QMLE}}},{\scriptscriptstyle{\textnormal{L}}}}=\sum_{a=1}^k c_a \widehat{\mu}_{n,a}^{{\scriptscriptstyle{\textnormal{QMLE}}},{\scriptscriptstyle{\textnormal{L}}}}$. The variance of $\widehat{\mu}_{n,c}^{{\scriptscriptstyle{\textnormal{QMLE}}},{\scriptscriptstyle{\textnormal{L}}}}$ (asymptotic variance of $\widehat{\mu}^{{\scriptscriptstyle{\textnormal{QMLE}}}}_{n,c}$) can be expressed as

No-harm-GR estimators

As discussed in freedman2008regression and cohen2020no, QMLE-GR estimators may perform worse than the baseline IPW estimator in terms of asymptotic variances under "misspecification" of the adjusting model. In this section, we consider GR estimators that do not increase asymptotic variances regardless of the configurations of potential outcomes and covariates. Recall that for a parameter of interest $\mu_{n,c}=c'\frac{1}{n}\mathbf{1}'y$, the asymptotic variance of a QMLE-GR estimator is

There is no guarantee that this asymptotic variance is smaller than the variance of the IPW estimator $ \frac{1}{n}c' \mathbf{1}'\operatorname{diag}(y)\hspace{1pt}\mathbf{\Omega}\hspace{1pt} \operatorname{diag}(y)\mathbf{1} c$, which is the asymptotic variance of the IPW estimator. To guarantee that our estimation strategy does not result in any harm, we can further define a multiplicative constant $\alpha_n$ that solves the problem

Note that the minimum is guaranteed to be no larger than $\frac{1}{n}c' \mathbf{1}'\operatorname{diag}(y)\hspace{1pt}\mathbf{\Omega}\hspace{1pt} \operatorname{diag}(y)\mathbf{1} c$ since $\alpha=0$ is in the choice set. The analytical expression for $\alpha_n^c$ is

for which we can construct a feasible consistent estimator:

Algorithm (ref) gives a formal template for constructing the No-harm GR estimators.\footnote{ This class of estimators is motivated by cohen2020no's no-harm estimator for a two-arm completely randomized design. Their estimator is different from ours because they use the outputs of QMLE imputations (i.e., $f^a(x_i,\widehat{\theta}_n)$) as regressors for an interacted linear regression model and apply lin2013agnostic's results. In general, one can design no-harm estimators that are more flexible than the No-harm-GR estimators considered here. For example, one can define the multiplicative constants separately for each arm or use the QMLE imputations of arms as regressors, as in cohen2020no. We will discuss one class of such estimators in the section below.}

Statistical guarantees are derived under the following additional assumptions.

Define $\theta_n$ as the target QMLE coefficient as before and $\alpha_n^c$ as in ((ref)). Define the infeasible estimator:

for $a\in[ k]$ and $\widehat{\mu}_{n,c}^{{\scriptscriptstyle{\textnormal{NOH}}},{\scriptscriptstyle{\textnormal{L}}}}=\sum_{a=1}^k c_a \widehat{\mu}_{n,a}^{{\scriptscriptstyle{\textnormal{NOH}}},{\scriptscriptstyle{\textnormal{L}}}}$. The variance of $\widehat{\mu}_{n,c}^{{\scriptscriptstyle{\textnormal{NOH}}},{\scriptscriptstyle{\textnormal{L}}}}$ (asymptotic variance of $\widehat{\mu}^{{\scriptscriptstyle{\textnormal{NOH}}}}_{n,c}$) can be expressed as

Denote the $l_1$-induced matrix norm by $\left\lvert\kern-0.25ex\left\lvert\kern-0.25ex\left\lvert A \right\rvert\kern-0.25ex\right\rvert\kern-0.25ex\right\rvert_1=\max_{j\in [n]}\{\sum_{i=1}^n |a_{ij}|\}$.

Optimal GR estimators (Opt-GR)

QMLE-GR estimators and No-harm-GR estimators both elicit $\theta_n$ from other criteria. Instead of requiring $\theta_n$ to be the minimizer of a pseudo-criterion, we can choose $\theta_n$ to directly minimize the implied asymptotic variances. That is, we define

The coefficient $\theta_n$ is optimal in the sense that the implied asymptotic variance $\frac{1}{n}c' \mathbf{1}'\operatorname{diag}\left(y-f(\theta)\right)\hspace{1pt}\mathbf{\Omega}\hspace{1pt} \operatorname{diag}\left(y-f(\theta)\right)\mathbf{1} c$ is the smallest when compared with estimators using the same class of parametric models for adjustments.\footnote{If we replace the first-order design matrix $\mathbf{\Omega}$ with an identified variance bound $\widetilde{\mathbf{\Omega}}$, the implied $\theta_n$ will instead minimize the asymptotic variance bound.} However, it is impossible to construct a sample analog to the criterion above because some entries in $y-f(\theta_n)$ are never observed simultaneously. This is the same problem we encountered in variance estimation, where some pairs of potential outcomes can never be observed simultaneously.

To identify and estimate $\theta_n$, we modify the criterion to allow for a feasible sample analog. Specifically, we define:

and its sample analog:

The criterion $ \mathcal{L}_n^{{\scriptscriptstyle{\textnormal{Opt}}}}(\theta)$ differs from the original criterion ((ref)) by the term $\frac{1}{n^2}c' \mathbf{1}'\operatorname{diag}(y)\hspace{1pt}\mathbf{\Omega}\hspace{1pt} \operatorname{diag}(y)\mathbf{1} c$. This term does not depend on the parameter value $\theta$ and hence minimizers for ((ref)) and ((ref)) are identical. Furthermore, unlike ((ref)), $\mathcal{L}_n(\theta)$ avoids the joint non-observability issue caused by the term $\frac{1}{n^2}c' \mathbf{1}'\operatorname{diag}(y)\hspace{1pt}\mathbf{\Omega}\hspace{1pt} \operatorname{diag}(y)\mathbf{1} c$ and hence admits a valid sample analog $\widehat{\mathcal{L}}^{{\scriptscriptstyle{\textnormal{Opt}}}}(\theta)$.

The Opt-GR estimators are constructed according to Algorithm (ref).

Theoretical guarantees can be derived under the standard identification assumptions, with additional standard regularity assumptions (Assumption (ref)) which we include in Appendix (ref). Though the identification assumption is standard, it requires careful consideration in our setting, so we highlight it here. For a set $C \subset \mathbb{R}^s$, we denote its boundary under the standard Euclidean topology by $\textrm{Bd}\left(C\right)$.

Define $\widehat{\mu}_{n,a}^{{\scriptscriptstyle{\textnormal{Opt}}},{\scriptscriptstyle{\textnormal{L}}}}$, $\widehat{\mu}_{n,c}^{{\scriptscriptstyle{\textnormal{Opt}}},{\scriptscriptstyle{\textnormal{L}}}}=\sum_{a=1}^k c_a\widehat{\mu}_{n,a}^{{\scriptscriptstyle{\textnormal{Opt}}},{\scriptscriptstyle{\textnormal{L}}}}$, and $ \text{\textnormal{Var}}(\widehat{\mu}_{n,c}^{{\scriptscriptstyle{\textnormal{Opt}}},{\scriptscriptstyle{\textnormal{L}}}})$ as in ((ref)) and ((ref)), but with $\theta_n$ replaced by $\theta_n\equiv\arg\min_{\theta\in\Theta} \mathcal{L}_n(\theta)$ where $\mathcal{L}_n(\theta)$ is defined in ((ref)).

We end this section by noting that one can combine the insights of No-harm GR and Opt-GR estimators to form a class of Optimal-Imputed (Opt-I) GR estimators. This class of estimators is Opt-GR estimators with a linear imputation function. Instead of using the full set of covariates, the linear imputation function uses a single covariate imputed by a QMLE model (for example, a flexible ML model).\footnote{We note that the coefficient estimator for the optimal-GR estimator with linear imputation functions admits a closed-form solution.} We find that this class of estimators has good finite-sample performance in our simulations. We detail the constructions below. Recall the definition of the intercept matrix $\mathbf{1}$ in equation ((ref)) and the stack of imputations $f(\theta)$ defined before Section (ref). The assumptions, asymptotic theory, and variance-bound estimation for the Opt-I estimators are left to Appendix (ref) due to space constraints.

Data Application

We now demonstrate how the results from earlier sections can be applied in practice. We study a network experiment based on the data in cai2015social. Section (ref) describes the background and the dataset. Section (ref) uses the comlexity metric $\sigma_{\max}\left(\mathbf{\Omega}^c\right)$ to understand the strengths and weaknesses of different designs in this setting. Section (ref) introduces simulation designs and discusses simulation results.

Background and Dataset

cai2015social examines how social networks influence weather insurance adoption in rural China. The unit of the experiment is a household and the primary outcome of interest is weather insurance adoption, a binary variable indicating whether a household purchases insurance after an information session. Some household characteristics are observed, including demographics, rice production, income, and past experiences with natural disasters. The social network information is collected through a friend-nomination survey.\footnote{For more detailed information, see Section II-B of cai2015social.}

Households were randomly assigned to two rounds of information sessions. In each round, a simple information session and an intensive information session were held simultaneously. The two rounds of information sessions were three days apart. Households were randomized into four treatment arms: first-round simple, first-round intensive, second-round simple, and second-round intensive. The network effects on the insurance take-ups are measured by the average differences in insurance purchase decisions among second-round participants with different numbers of friends who were invited to the first round.

The experimental design is a village-level stratified randomization. Households were stratified according to household sizes and areas of rice production per capita. The computer code used for randomization in cai2015social is not immediately available.\footnote{A brief description of the randomization procedure is given in footnote 8 of cai2015social.} As a result, we decide to construct the randomization procedures based on the description of the paper and the details can be found in Appendix (ref). Any implications we draw below should not be related to the original paper.

We briefly describe the social networks used in the data. The friend nomination graph is a directed graph with edges pointing from the nominators to the nominees. The network has 4806 households and 41 disjoint weak components.The average outdegree is 3.5, the maximum indegree is 18, and the average within-village path length is 2.7. We plot the second largest component of the network in Figure (ref). The network has clear community structures, as most friendship ties form within natural villages. {

} There are 12 exposure mappings considered in the paper. We include a subset of them in Table (ref). Exposures 8–12 are defined analogously to exposures 3–7, replacing the SRS with the SRI. Motivated by Column (5) in Table 2 of cai2015social, we focus on the comparisons between exposures 3v4, 3v5, 3v6, 4v5 and 4v6, which we treat as the key parameters of interest to measure the effects of social networks.

Use $\sigma_{\max}\left(\Omega^c\right)$ to understand different designs

In this section, we compare three designs using our proposed measure $\sigma_{\max}\left(\Omega^c\right)$ from Section (ref). The purpose of this demonstration is to show how to use the measure to understand the relative strengths and weaknesses of various designs. We consider multiple experimental designs: A) a finely-stratified village-level randomization, B) a village-level randomization, and C) Bernoulli designs with varying assignment probabilities.

The finely-stratified village-level randomization (Design A) is inspired by the experimental design in cai2015social. We partition households in each natural village into four groups based on their household sizes and rice production areas. Unis within each stratum are completely randomized into four treatment arms. For details, see Appendix (ref). For the village-level randomization (Design B), we randomly assign households within each village to four arms with proportions (0.1,0.1,0.4,0.4). For Bernoulli designs, we consider different assignment probabilities to four treatment arms: (1/4,1/4,1/4,1/4), (1/5,1/5,3/10,3/10), (1/6,1/6,4/12,4/12) and (1/9,2/9,4/12,4/12). We label them Designs C.1, C.2, C.3 and C.4, respectively.

Table (ref) below documents $\sigma_{\max}\left(\Omega^c\right)$ values for different designs. We make two comments about the tables:

Simulation Design

We report simulation results for comparing exposures 3 vs.4 and comparing exposures 3 vs.6 using Design C.2 and Design C.4. The goal of this simulation is to illustrate implications of the complexity metric $\sigma_{\max}\left(\mathbf{\Omega}^c\right)$ as well as the behavior of the various estimators.

We compare the following 7 estimators: 1) IPW estimator, 2) inverse-probability weighted WLS estimator (WLS), 3) QMLE-GR estimator with a logit model (Logit),\footnote{We set $\omega_{ai}=\pi_{ai}$ where $\omega_{ai}$ is defined in equation ((ref)). This is to mimic the exercise where researchers estimate a logit model without any weighting.}, 4) Opt-GR estimator with a linear model (Opt Linear), 5) Opt-GR estimator with a logit model (Opt Logit), 6) Opt-I GR estimator with imputations using the WLS model (Opt-I WLS), and 7) Opt-I GR estimator with imputations using the logit model (Opt-I Logit). All adjustment models have a separate intercept for each arm and the same coefficients on the covariates.

For each comparison of two exposures, we impute the potential outcomes in two ways. In the first simulation scenario, we impute the potential outcomes using a logistic model with coefficients estimated from the data. In this scenario, barring finite-sample issues, QMLE-GR, Opt-GR, and Opt-I GR estimators are expected to work similarly well and show improvement over the baseline IPW estimator. In the second simulation scenario, we impute the potential outcomes such that the Opt-GR estimators will have efficiency gains over the QMLE-GR estimators. This scenario is used to demonstrate the theoretical guarantee for the Opt-GR and Opt-I GR estimators. We refer to the first simulation scenario as Sim-Impute and the second simulation scenario as Sim-Optimal. More details for data constructions, imputations, and implementations can be found in Appendix (ref).

Simulation Results

Formal and complete simulation results are included in Appendix (ref). We high a subset of the results in Table (ref) and Table (ref) for discussion:

Additional Results in Section (ref)

We define a few more estimators.

The following lemma gives conditions under which the CI and MI estimators are consistent estimators of the average potential outcomes. These results are stated in an unpublished work middleton2021unifying3. See also brewer1979class and wright1983finite.

The following lemma gives conditions under which a WLS estimator is algebraically equivalent to a GR estimator:

Additional Results in Section (ref)

This section contains additional results for the Opt-I estimators. Define $X=

\in\mathbb{R}^{kn\times(k+1)}$ and $\mathbf{\Omega}^c= \operatorname{diag}(c'\mathbf{1})\mathbf{\Omega} \operatorname{diag}(c'\mathbf{1})\in \mathbb{R}^{kn\times kn}$.

We define:

for $a\in[k]$ and $\theta_n\equiv\arg\min_{\theta\in\Theta} \mathcal{L}_n(\theta)$ defined in ((ref)). Define $\hat{\mu}_{n,c}^{{\scriptscriptstyle{\textnormal{Opt-I}}},{\scriptscriptstyle{\textnormal{L}}}}=\sum_{a=1}^k c_a\hat{\mu}_{n,a}^{{\scriptscriptstyle{\textnormal{Opt-I}}},{\scriptscriptstyle{\textnormal{L}}}}$. The variance of $\hat{\mu}_{n,c}^{{\scriptscriptstyle{\textnormal{Opt-I}}},{\scriptscriptstyle{\textnormal{L}}}}$ can be expressed as

Additional Results for the Network Experiment in Section (ref)

We provide lower-level conditions for the network experiments presented in Section (ref). These conditions lead to the properties $\sigma_{\max}\left(\mathbf{\Omega}\right)=O(1)$, $\sigma_{\max}((\widetilde{\mathbf{\Omega}}\otimes \widetilde{\mathbf{\Omega}})\circ \mathbf{S})=o(n)$ and $\sigma_{\max}\left(\widetilde{ \mathbf{\Omega} }_{\hspace{-.6mm}{}{/}}{}_{\scriptscriptstyle \hspace{-.6mm}\mathbf{p}}\right)=O(1)$, which are required by Corollary (ref) for $\sqrt{n}$-estimation and Corollary (ref) for consistent variance bound estimation.

One can verify that the Assumptions (ref) and (ref) are satisfied for the experimental design considered in Example (ref) when the network has a bounded degree uniformly in $n$.\footnote{ In some social network settings, there may exist nodes of very large degrees. This may cause two problems: 1) the assignment probability becomes too small, and/or 2) the dependence between exposure maps becomes too strong. In these cases, one may wish to restrict the parameters of interest to subgroups (e.g., people with low degrees) and consider specific exposure mappings for precise estimates (e.g., the sample average treatment effects condition on the event that all the high-degree nodes are assigned to treatment). }

These two assumptions are standard in the literature. Under the two assumptions, both the first-order design matrix $\mathbf{\Omega}$ and the AS bound matrix $\tilde{\mathbf{\Omega}}^{{\scriptscriptstyle{\textnormal{AS}}}}$ have bounded $l_1$-induced norms, because each column contains only finitely many nonzero entries with bounded magnitudes uniformly in $n$. The fact that $\sigma_{\max}((\widetilde{\mathbf{\Omega}}\otimes \widetilde{\mathbf{\Omega}})\circ \mathbf{S})=o(n)$ and $\sigma_{\max}\left(\widetilde{ \mathbf{\Omega} }_{\hspace{-.6mm}{}{/}}{}_{\scriptscriptstyle \hspace{-.6mm}\mathbf{p}}\right)=O(1)$ are also satisfied by Proposition 6.2 in aronow2017estimating and by Lemma (ref).

Recall that $p$ is the dimension of the pretreatment covariates. We consider the case of a same-slope adjustment and consider two types of imputation functions: linear models $f^a(x_i,\theta)=\gamma^a+x_i'\beta$, $a\in[k]$ and logistic models $f^a(x_i,\theta)=\frac{\exp(\gamma^a+x_i'\beta)}{1+\exp(\gamma^a+x_i'\beta)}$, $a\in [k]$. For the QMLE-GR estimator with the linear model, the finite population criterion and the sample equivalent are

where $\theta=(\{\gamma^a\}_{a=1}^k,\beta)\in\mathbb{R}^{k+p}$. For the QMLE-GR estimator with the logistic regression model, they are

and

where $\theta=(\{\gamma^a\}_{a=1}^k,\beta)\in\mathbb{R}^{k+p}$. Note we can also define

and

For the Opt-GR estimators, we form the criterion using the setup in Section (ref) with the first-order design matrix $\mathbf{\Omega}$. We choose the AS bound for a bounding matrix. The plug-in variance bound estimators are constructed using the formulae in Theorems (ref), (ref), and (ref). For the Opt-GR estimator with the logistic model, we need more moment assumptions in order to satisfy Assumption (ref)-(i)-(b).

The following theorems specialize Theorems (ref), (ref), and (ref) in the network experiment setting.

Proofs for this theorem can be found in Appendix (ref).

Auxiliary Lemmas

This section proves several auxiliary lemmas instrumental in the proofs below. Let $\|\cdot\|_2$ denote the Frobenius norm if applied to a matrix and the $l_2$ vector norm if applied to a vector.

Let $z\in\mathbb{R}^{kn}$ be an arbitrary vector. We define the IPW estimator:

The following lemma shows that, in order to establish convergence for a symmetric matrix estimator (of a fixed dimension), it suffices to consider its bilinear forms.

We now state a tensor inequality. Consider a fourth-order $n$-dimensional tensor. We denote it as $\mathbf{A}=[a_{ijkl}]\in\mathbb{R}^{n\times n\times n \times n}$. We shall understand it as a multi-linear function:

where for $x,y,z,a\in\mathbb{R}^n$,

Consider the following maximization problem:

We denote its optimal value as $\sigma_{\max}(\mathbf{A})$. This optimal value exists because we are optimizing a continuous function over a compact set. Note further the problem

has the same solution as the problem above because we can always take the negative of one of the vectors.

Note for any vectors $w,x,y,z\in\mathbb{R}^n$, we have:

The following lemma bounds $\sigma_{\max}(\mathbf{A})$. We define quantities:

and analogously,

and similarly for $\|A\|_{-k}$ and $\|A\|_{-l}$. We define:

Denote $\|z\|_4^4=\sum_{i,j}z^4_{ij}$, $\|z\|_1=\sum_{i,j}|z_{ij}|$ and $\|z\|_{\infty}=\max_{i,j}|z_{ij}|$ for a matrix $z$.

Let $\left\lvert\kern-0.25ex\left\lvert\kern-0.25ex\left\lvert \mathbf{A} \right\rvert\kern-0.25ex\right\rvert\kern-0.25ex\right\rvert_1$ denote the $l_1$-induced matrix norm, where $\left\lvert\kern-0.25ex\left\lvert\kern-0.25ex\left\lvert \mathbf{A} \right\rvert\kern-0.25ex\right\rvert\kern-0.25ex\right\rvert_1=\max_{j\in[n]}\{\sum_{i=1}^n |a_{ij}|\} $, for an arbitrary $\mathbf{A}\in\mathbb{R}^{m\times n}$.

Let $\Theta$ be a compact set in a finite dimensional Euclidean space, and $B(\theta,\delta)$ denote a closed ball in $\Theta$ of radius $\delta \geq 0$ (with the $l_2$ norm) centered at $\theta$. The following lemma adapts Theorem 1 andrews1992generic to our setting.

The following lemma is the standard consistency proof for GMM estimators.

The following lemma uses notations in Section (ref).

Proofs for the results in Section (ref)

Proof of Lemma (ref)

Notice that if $\frac{1}{n}\sum_{i=1}^n x_{si}=0$ for all $s\in [p]$, we have $ \frac{1}{n}\mathbf{1}'X =

$. Hence $

b^{{\scriptscriptstyle{WLS}}}=\mathbf{1}'Xb^{{\scriptscriptstyle{WLS}}}$. Hence $\frac{1}{n}\mathbf{1}'\left(y - Xb^{{\scriptscriptstyle{WLS}}}\right)=0$ if columns of $\mathbf{1}$ are in the column space of $\boldsymbol{\pi}\OmX$ since

for some $t\in\mathbb{R}^{(k+p)\times k}$ such that $\mathbf{1}=\pi\mathbf{\omega} X t$, by the first-order conditions of a WLS estimator.

Proof of Theorem (ref)

Proof of Corollary (ref)

Proof of Theorem (ref)

Proof of Corollary (ref)

Proof of Theorem (ref)

Proof of Corollary (ref)

We prove only for the GR estimators. The corollary can be proved in an analogous way for other estimators. The plug-in variance bound estimator for a GR estimator is:

Using the notation in Theorem (ref), we identify $z=\mathbf{1}'\operatorname{diag}(y-X b^{{\scriptscriptstyle{\textnormal{WLS}}}}_n)\in\mathbb{R}^{kn\times k}$ and $\widehat{z}={D} \mathbf{1}'\operatorname{diag}(y-X \widehat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_n)\in\mathbb{R}^{kn\times k}$. Note $X(\widehat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_n-b^{{\scriptscriptstyle{\textnormal{WLS}}}}_n)\in\mathbb{R}^{kn}$. We have

where $\frac{1}{n}\|X\|_2^2=O(1)$ by Assumption (ref) and $\|\widehat{b}^{{\scriptscriptstyle{\textnormal{WLS}}}}_n-b^{{\scriptscriptstyle{\textnormal{WLS}}}}_n\|_2^2=O_p\left(\sqrt{\sigma_{\max}\left(\mathbf{\Omega}\right)}n^{-\frac{1}{2}}\right)$ by Lemma (ref). $\frac{1}{n}\|z\|_2^2=O(1)$ by Assumption (ref) and Lemma (ref).

Proof of Theorem (ref)

Based on the premises, we have $\widehat{\nu}_n-\widehat{\nu}_n^{{\scriptscriptstyle{\textnormal{L}}}}=o_p\left(\sqrt{\sigma_{\max}\left(\mathbf{\Omega}\right)}/n^{0.5}\right)$. Because $\text{\textnormal{Var}}\left(c'\widetilde{\nu}_n^{{\scriptscriptstyle{\textnormal{L}}}}\right)\succeq \sigma_{\max}\left(\mathbf{\Omega}\right)/n$, we have $\left(c'\widehat{\nu}_n-c’\widehat{\nu}^{{\scriptscriptstyle{\textnormal{L}}}}_n\right)/\sqrt{\text{\textnormal{Var}}\left(c'\widehat{\nu}_n^{{\scriptscriptstyle{\textnormal{L}}}}\right)}=o_p(1)$.

We have

Hence:

for any positive $\epsilon$. The second term converges to 0. Since $q(\alpha)$ is continuous, for sufficiently small $\epsilon$, there exists a $\delta$ such that $q(\alpha-\delta)\leq (1+\epsilon)^{-1}q(\alpha)$. Hence we have,

Since the $\epsilon$ can be made arbitrarily small, the $\delta$ can be made arbitrarily small as well, and we have

Additional Assumptions

We use $\nabla$ to denote the total differentiation operator. For a function $f:\mathcal{X}\subset \mathbb{R}^k\to\mathbb{R}$, $\nabla_x f$ denotes the gradient function of $f$ (if it exists), $\nabla_{xx'}f$ denotes the Hessian function of $f$ and so on. We use the partial derivative notation $\frac{\partial}{\partial x}f$ to denote the partial derivative $f$ with respect to a particular argument $x$. For a set in $\Theta\subset \mathbb{R}^s$, we use the notation $\operatorname{Bd}(\Theta)$ to denote its boundary with respect to the standard topology of a Euclidean space.

We make the following assumptions to study the probabilistic properties of QMLE-GR estimators. Let $\Theta$ be a set in a finite-dimensional Euclidean space. We define the distance from a point $\theta$ to a set $\Theta$ as $d(\theta,\Theta)=\inf_{\tilde{\theta}\in \Theta}\|\tilde{\theta}-\theta\|_2$. The following assumptions are used to obtain consistency and $\sqrt{n}$ convergence of $\hat{\theta}_n$. For a vector/matrix/tensor $A$, we use $\|A\|_1$ to denote the sum of the absolute values of the entries in $A$.

The following assumptions on $f^a(\cdot,\theta)$, $a\in [k]$ are used to obtain the $\sqrt{n}$ equivalence of the QMLE-GR estimators to their asymptotic linear expansions.

Proofs for the results in Section (ref)

Proof of Theorem (ref)

We first check the stochastic equicontinuity for large $n$. Notice by Assumption (ref)-(iv), $\sigma_{\max}\left(\mathbf{\Omega}\right)=O(1)$ and Lemma (ref) we have pointwise convergence for the criterion function: for each $\theta\in\Theta$,

By Assumption (ref)-(v), $\mathcal{L}_n(\theta)$ is continuous uniformly over $\theta\in\Theta$ and for large $n$. Then,

First note that the second term is a degenerate probability event: it happens with probability 0 for sufficiently small $\delta$. For the first term, we take a $2\delta$-covering $\mathcal{N}_{2\delta}$ of $\Theta$. Note $\mathcal{N}_{2\delta}$ is a finite set by Assumption (ref)-(i). By the triangle inequality, for each $\theta$ there exists a $\tilde{\theta}\in\mathcal{N}_{2\delta}$ such that $d(\theta',\tilde{\theta})<2\delta$ for all $\theta'\in B(\theta,\delta)$. Thus we can bound the first term:

where the last line is by the Markov inequality and Assumption (ref)-(ix). By setting $\delta$ small enough, we prove the desired inequality with Assumption (ref)-(v). With Lemma (ref) and Lemma (ref), we conclude $\widehat{\theta}_n-\theta_n=o_p(1)$. We now show $\widehat{\theta}_n-\theta_n=o_p(n^{-\frac{1}{2}})$.

By Assumption (ref)-(ii), $\widehat{\theta}_n$ is in the interior of $\Theta$ with probability approaching one. Hence with probability approaching one, $\widehat{\theta}_n$ satisfies:

A Taylor expansion around $\theta=\theta_n$ yields:

The last term is justified by bounding the higher-order reminder terms as follows: take the $t$th entry of the gradient $ \nabla_{\theta}\widehat{\mathcal{L}}_n$. This row corresponds to the partial derivative of $\widehat{\mathcal{L}}_n$ with respect to the $t$th parameter $\theta_t$. Its Taylor expansion has the form:

where $\tilde{\theta}_n$ is between $\theta_{n}$ and $\hat{\theta}_{n}$. As $\widehat{\theta}_n$ enters $B(\theta_n,\epsilon)$ with probability one and by the subadditivity of the spectral norm $\left\lvert\kern-0.25ex\left\lvert\kern-0.25ex\left\lvert \cdot \right\rvert\kern-0.25ex\right\rvert\kern-0.25ex\right\rvert_2$, we have,

where for ((ref)) we use the fact that the spectral norm is dominated by the Frobenius norm, and for ((ref)) we used the the equivalence of the $l_1$ and $l_2$ norm, Markov inequality, Assumption (ref)-(vii), Assumption (ref)-(viii) and $\sigma_{\max}\left(\mathbf{\Omega}\right)=O(1)$. Also by Assumption (ref)-(vi), Assumption (ref)-(viii) and Lemma (ref), $\frac{1}{n}\sum_{a,i} \omega_{ai}\frac{{D}_{ai}}{\boldsymbol{\pi}_{ai}}\nabla_{\theta\theta}g^a( y_{i}(a),x_i,\theta_n)$ converges in probability to $\frac{1}{n}\sum_{a,i}\omega_{ai} \nabla_{\theta\theta}g^a( y_{i}(a),x_i,\theta_n)$. Hence it is invertible with probability approaching one by Assumption (ref)-(ix). Using the above argument and the fact that $\widehat{\theta}_n-\theta_n=o_p(1)$, we rearrange ((ref)) and ((ref)) to get to

where for the last line we used the first order condition $\frac{1}{n}\sum_{i=1}^n\omega_{ai}\nabla_\theta g^a( y_{i}(a),x_i,\theta_n)=0$, Assumption (ref)-(vi), Assumption (ref)-(viii), Lemma (ref) and $\sigma_{\max}\left(\mathbf{\Omega}\right)=O(1)$.

Because $\sigma_{\max}((\widetilde{\mathbf{\Omega}}\otimes \widetilde{\mathbf{\Omega}})\circ \mathbf{S})=o(n)$, $\sigma_{\max}\left(\widetilde{ \mathbf{\Omega} }_{\hspace{-.6mm}{}{/}}{}_{\scriptscriptstyle \hspace{-.6mm}\mathbf{p}}\right)=O(1)$ and $n\text{\textnormal{Var}}(\widehat{\mu}_{n,c}^{{\scriptscriptstyle{\textnormal{QMLE}}},{\scriptscriptstyle{\textnormal{L}}}})\geq c_{\ref{Thm:QMLE}}$ uniformly for all large $n$, by Assumption (ref), Lemma (ref), Lemma (ref) and Lemma (ref), and the fact that

where $f(\theta)$ is the vector of imputed outcomes of all arms using coefficient $\theta$, we have,

and

The remainder of the proofs is similar to that for Theorem (ref).

Proof of Theorem (ref)

We first show that

The numerator is upper bounded by

is bounded above uniformly in $n$ by Assumption (ref), $\sigma_{\max}\left(\mathbf{\Omega}\right)=O(1)$, and Assumption (ref)-(ii). The denominator is bounded below uniformly in $n$ by Assumption (ref). Thus the imputation functions $\alpha f^a(,\theta)$, $a\in [k]$, satisfy Assumption (ref).\\ Estimator for the denominator converges at a $\sqrt{n}-$rate by Lemma (ref),

by Lemma (ref). The numerator also converges at a $\sqrt{n}-$rate. We first show:

Denote $c'\mathbf{1}' f(\theta_n)$ by $f^c(\theta_n)\in\mathbb{R}^{kn}$ and $c'\mathbf{1}\operatorname{diag}(y)$ by $y^c\in\mathbb{R}^{kn}$. The term $\frac{1}{n}c'\mathbf{1}'f(\theta_n)\mathbf{\Omega} \boldsymbol{\pi}^{-1} ({D}-\boldsymbol{\pi})\operatorname{diag}(y)\mathbf{1} c$ is $O_p(n^{-\frac{1}{2}})$ by noticing:

and by Lemma (ref), Assumption (ref), and Assumption (ref). Together with Assumption (ref) and the fact that $\left\lvert\kern-0.25ex\left\lvert\kern-0.25ex\left\lvert \mathbf{\Omega} \right\rvert\kern-0.25ex\right\rvert\kern-0.25ex\right\rvert_1=O(1)$, this implies $\hat{\alpha}_n^c-\alpha_n^c=O_p(n^{-\frac{1}{2}})$. The remaining proof is the same as the proof in Theorem (ref).

Proof of Theorem (ref)

Notice first:

We only need to check the $\sqrt{n}$-consistency $\widehat{\theta}_n$. The rest proofs are identical to those of Theorem (ref). We first check the pointwise convergence of the criterion, namely

pointwise in $\theta\in\Theta$. We adopt the notation $y^c$ and $f^c\left(\theta\right)$ in the proof of Theorem (ref).

First note that $ \widehat{\mathcal{L}}_n(\theta)$ is an unbiased estimator of $ \mathcal{L}_n(\theta)$. The variance of $ \widehat{\mathcal{L}}_n(\theta)$ is:

where for the last inequality we use Lemma (ref). Because $\left\lvert\kern-0.25ex\left\lvert\kern-0.25ex\left\lvert \mathbf{\Omega} \right\rvert\kern-0.25ex\right\rvert\kern-0.25ex\right\rvert_1=O(1)$ and $\sigma_{\max}\left(\mathbf{\Omega}\right)=O(1)$ and under Assumption (ref) and Assumption (ref)-(i)-(a), the term above is $o_p(1)$. \\ Next we check stochastic equicontinuity as in the proof of Theorem (ref):

for a constant $C$ bounded above uniformly for large $n$ by Assumption (ref) and Assumption (ref)-(ii)-(b). Using the same steps as in Theorem (ref), we can prove stochastic equicontinuity upon noticing

as $\|\theta_1-\theta_2\|_2\to 0$. Further with Assumption (ref)-(ii) and by Lemma (ref) and Lemma (ref), we have $\widehat{\theta}_n-\theta_n=o_p(1)$.

Now we establish the rate of convergence. We define $\frac{\partial}{\partial \theta_t}f^c(\theta)\equiv c'\mathbf{1}'\operatorname{diag}(\frac{\partial}{\partial \theta_t} f(\theta))\in \mathbb{R}^{kn}$. Let $\nabla_{\theta} f^c(\theta)\in \mathbb{R}^{kn\times s}$ denote the column stacks of $\frac{\partial}{\partial \theta_t}f^c(\theta), t\in [s]$ . By Assumption (ref)-(ii), with probably approaching one, $\widehat{\theta}_n$ is in the interior of $\Theta$ and satisfies the first order condition:\footnote{We omit the constant factor $2$.}

For the $t$th entry of the equalities above, it can be written as:

We expand:

We study (A), (B), (C), and (D). For (A), we have the decomposition:

where the last line is by Assumption (ref) , Assumption (ref)-(iii)-(b) and that $\sigma_{\max}\left(\mathbf{\Omega}\right)=O(1)$.

To study (B), we have the decomposition:

by Assumption (ref), Assumption (ref)-(i)-(a), Assumption (ref)-(iii)-(b) and $\sigma_{\max}\left(\mathbf{\Omega}\right)=O(1)$.

To study (C), we have the decomposition,

by Assumption (ref)-(i)-(b) and Assumption (ref)-(iii)-(a) and $\sigma_{\max}\left(\mathbf{\Omega}\right)=O(1)$.

For $(D)$, we have:

by Assumption (ref)-(ii)-(b), Assumption (ref)-(ii)-(c) and $\sigma_{\max}\left(\mathbf{\Omega}\right)=O(1)$.

To summarize, for the $t$th column of the equalities in equation ((ref)), we have:

Define $\widehat{\mathcal{H}}_n$ to be the column stack of the vectors $\frac{1}{n}\left(y^c{D}\boldsymbol{\pi}^{-1}-f^c\left(\theta_n\right)\right)\mathbf{\Omega} \nabla_{\theta}\frac{\partial}{\partial \theta_t}f^c\left(\theta_n\right)$, $t\in [s]$.

Stacking the previous expressions together, we have the expression:

Define,

where $\mathcal{H}_n$ is defined in Assumption (ref)-(iv). Note that $\text{\textnormal{E}}\left[\widehat{b}_n\right]=0$ by Assumption (ref)-(ii), and $\widehat{b}_n=O_p(n^{-\frac{1}{2}})$ by Assumption (ref), Assumption (ref)-(i)-(a), Lemma (ref) and the fact that $\left\lvert\kern-0.25ex\left\lvert\kern-0.25ex\left\lvert \mathbf{\Omega} \right\rvert\kern-0.25ex\right\rvert\kern-0.25ex\right\rvert_1=O(1)$.

We have $\widehat{A}_n-A_n=o_p(1)$ by Assumption (ref), Assumption (ref)-(i), Lemma (ref) and the fact that $\left\lvert\kern-0.25ex\left\lvert\kern-0.25ex\left\lvert \mathbf{\Omega} \right\rvert\kern-0.25ex\right\rvert\kern-0.25ex\right\rvert_1=O(1)$. Further, we have $A_n=O(1)$ by Assumption (ref)-(i), and $A_n^{-1}=O(1)$ and $\widehat{A}_n^{-1}=O_p(1)$ by Assumption (ref)-(iv). Hence we have,

Rest of the proofs are similar to those in Theorem (ref).

Proof of Theorem (ref)

The proof is identical to that of the Theorem (ref).