\section{Motivation: randomization-based causal inference}
\label{sec:motivation}

\subsection{Existing results}

In a seminal paper, \citet{neyman1923application} introduced the notation of potential outcomes to define causal effects.  More importantly, he also proposed a framework for statistical inference of causal effects based on the design of the experiment. In particular, \citet{neyman1923application} considered an experiment with $N$ units and $Q$ treatment arms, where the number of units under treatment $q$ equals $N_q$, with $\sum_{q=1}^Q N_q= N$. Corresponding to treatment level $q $, unit $i$ has the potential outcome $Y_i(q)$, where $i=1, \ldots, N$ and $q = 1, \ldots, Q$. Despite its simplicity, the following completely randomized experiment has been widely used in practice and has generated rich theoretical results. Definition \ref{def:cr} below characterizes the joint distribution of ${{\boldsymbol{Z}}} = (Z_1, \ldots, Z_N)$, where $Z_i \in \{1, \ldots, Q\}$ is the treatment indicator for unit $i$.

\begin{definition}[Complete randomization]\label{def:cr}
$ {\mathbb{P}}( {\boldsymbol{Z}} = {\boldsymbol{z}}) = N_1 ! \cdots N_Q! / N! $ for all ${\boldsymbol{z}} = (z_1, \ldots, z_N ) $ with $(N_1, \ldots, N_Q)$ units receiving treatment level $(1,\ldots,Q)$, respectively.
\end{definition}

\citet{neyman1923application} formulated complete randomization based on an urn model, which is equivalent to Definition \ref{def:cr}.  Under complete randomization, ${\boldsymbol{Z}}$ is from a random permutation of $N_1$ $1$'s, $\ldots $, $N_Q$ $Q$'s. The experiment reveals one of the potential outcomes, which is the observed outcome $  Y_i = Y_i(Z_i) = \sum_{q=1}^Q Y_i(q) \ind{Z_i = q}$ for each unit $i$.

In \citet{neyman1923application} 's framework, all potential outcomes are fixed and only the treatment indicators are random according to Definition \ref{def:cr}. \citet[][Chapter 9]{scheffe1959analysis} called it the {\it randomization model}. Under this model, it is conventional to call the resulting inference as {\it randomization inference} or {\it design-based inference}. It has become increasingly popular in both theory and practice \citep[e.g.,][]{kempthorne52, copas1973randomization, robins1988confidence, rosenbaum2002observational, hinkelmann07, freedman2008Aregression, freedman2008Bregression, lin2013agnostic, dasgupta2015causal, imbens15, ATHEY201773, fogarty2018regression}.

A central goal in  \citet{neyman1923application} 's framework is to use the observed data $(Z_i, Y_i)_{i=1}^N$ to make inference of causal effects defined by the potential outcomes.  Define
$$
\overline{Y}(q) = N^{-1}\sum_{i=1}^N Y_i(q) ,\quad
S(q,q') =  (N-1)^{-1} \sum_{i=1}^N(Y_i(q)-\overline{Y}(q))(Y_i(q')-\overline{Y}(q'))
$$
as the average value of the potential outcomes under treatment $q$ and the covariance of the potential outcomes under treatments $q$ and $q'$, respectively.
Define the average potential outcome vector as $\overline{Y} = (\overline{Y}(1), \ldots,  \overline{Y}(Q))^\top  \in {\mathbb{R}}^Q$, and define the covariance matrix of the potential outcomes as $S = (S(q,q'))_{q,q'=1,\ldots, Q}$. The parameter of interest is a linear transformation of $\overline{Y}$:
\begin{align*}
    \gamma = F^\top \overline{Y}
\end{align*}
for a pre-specified $F = (f_{qh})\in{\mathbb{R}}^{Q\times H}$, which is called the contrast matrix when its columns are orthogonal to $(1,\ldots,1)^\top$. Despite the simple form of $\gamma$, it can answer questions from a wide ranges of applications.
For instance, \citet{neyman1923application} considered pairwise differences in means, and \citet{dasgupta2015causal} and \citet{mukerjee2018using} considered linear combinations of the mean vector. Recently, \citet{li2017general} unified the literature by studying the properties of the moment estimator for $\gamma$ under complete randomization. In particular, define the sample mean and variance of the observed  $Y_i(q)$'s as
\begin{align}\label{eqn:hS}
 {\hat{Y}}_q = N_q^{-1}\sum_{Z_i = q}Y_i,\quad {\hat{S}}(q,q) =  (N_q-1)^{-1} \sum_{Z_i = q} (Y_i - {{\hat{Y}}_q} )^2,
\end{align}
respectively. Define
\begin{align}\label{eqn:hY-hV}
{\hat{Y}} = ({\hat{Y}}_1,\cdots,{\hat{Y}}_Q)^\top\in{\mathbb{R}}^Q,\quad  {\hat{V}}_{\hat{Y}} = \diag{N_q^{-1}{\hat{S}}(q,q)}_{q\in[Q]} \in {\mathbb{R}}^{Q\times Q}
\end{align}
as the vector of sample averages and the diagonal matrix of the sample variances across all  arms, respectively. \citet{li2017general} showed that  $ {\hat{Y}}_q $ has mean and covariance
\begin{align}\label{eqn:mean-var-Gamma}
{\mathbb{E}}\{{\hat{Y}}_q \} = \overline{Y}  ,\quad
\Var{{\hat{Y}}_q} =  V_{\hat{Y}} = \diag{N_q^{-1}S(q,q)}_{q\in[Q]} - N^{-1}S,
\end{align}
and moreover, ${\hat{V}}_{\hat{Y}}$ is a conservative estimator for  $V_{\hat{Y}}$ in the sense that $ {\mathbb{E}}\{{\hat{V}}_{\hat{\gamma}}\} - V_{\hat{\gamma}}$ positive semi-definite. An immediate consequence is that
    \begin{align}\label{eqn:estimates}
        {\hat{\gamma}} = F^\top {\hat{Y}}, \quad {\hat{V}}_{\hat{\gamma}} = F^\top {\hat{V}}_{\hat{Y}} F
    \end{align}
are an unbiased point estimator for $    \gamma$ and a conservative covariance estimator for $V_{\hat{\gamma}} = \Var{{\hat{\gamma}}} =  F V_{\hat{Y}} F^\top$, respectively. \citet{li2017general} also used the established combinatorial or rank central limit theorems (CLTs) \citep{hajek1960limiting, hoeffding1951combinatorial, fraser1956vector} to prove the asymptotic Normality of $   {\hat{\gamma}} $ and the validity of the associated large-sample Wald-type inference, under certain regularity conditions.

\subsection{Open questions}
\label{sec::open-questions}

Despite the long history of \citet{neyman1923application}'s randomization model, the theory for randomization-based causal inference is far from complete. Technically, \citet{li2017general}'s review only covered the first regime \ref{regime:R1} below, and even there, finer results such as Berry--Esseen bounds (BEBs) have not been rigorously established for the most general setting (see \citet{wang2021rerandomization} for some recent results of BEBs for treatment-control experiments). For other regimes below, many basic results are missing in the literature. Table \ref{tab:many-settings} provides a list of important regimes and reviews the established and missing theoretical results. The discussion below gives more details.

\begin{table}[t]
\centering
\caption{Theoretical results for multi-armed experiments under the randomization model.
\ref{regime:R1}--\ref{regime:R4} correspond to nearly uniform designs by Definition \ref{def:uniform-design}, whereas  \ref{regime:R5} corresponds to non-uniform designs by Definition \ref{def:non-uniform-design}.
}
\label{tab:many-settings}
\begin{tabular}{M{1.0cm}M{1.0cm}M{3cm}M{6cm}}
\toprule
 Regime & $Q$   & $N_q$                    & CLT, variance estimation, and BEB              \\ \hline
   \ref{regime:R1} &  Small &   Large & CLT and variance estimation; no BEB  \\
\ref{regime:R2}        & Large & Large                    & Seems similar to \ref{regime:R1} but not studied \\
 \ref{regime:R3}        & Large & Small but $N_q\ge$ 2 & Not studied    \\
   \ref{regime:R4} &  Large &  $N_q = 1$ & Not studied; variance estimation is nontrivial \\ \hline
   \ref{regime:R5} &   \multicolumn{2}{c}{{Mixture of the above}} &  Not studied \\ \bottomrule
\end{tabular}%
\end{table}

\begin{enumerate}[label={({R\arabic*})}]
    \item\label{regime:R1} Small $Q$ and large $N_q$'s.  In this regime, the number of arms is small and the sample size in each arm is large. Asymptotically, as $N \to \infty$, we have that $Q$ is a fixed integer and $N_q/N\to e_q\in(0,1)$ for all $q=1,\ldots, Q$.
    \cite{li2017general} showed that, under \ref{regime:R1} and some regularity conditions on the potential outcomes, we have
\begin{equation}\label{eq::basic-asymptotic-inference}
        V_{\hat{\gamma}}^{-1/2}({\hat{\gamma}} - \gamma) \rightsquigarrow {\mathcal{N}}(0, I_H), \quad  N{\hat{V}}_{\hat{\gamma}} - N{\mathbb{E}}\{{\hat{V}}_{\hat{\gamma}}\} = o_{\mathbb{P}}(1),
\end{equation}
which ensure that the large sample Wald-type inference based on the Normal approximation is conservative.    \cite{li2017general}'s results are asymptotic. An important theoretical question is to quantify the finite-sample properties of ${\hat{\gamma}}$ by deriving non-asymptotic results.

    \item\label{regime:R2} Large $Q$ and large $N_q$'s. In this regime, each arm has adequate units for the variance estimation, but the number of arms is also large.  Asymptotically, as $N \to \infty$, we have $Q\to \infty$ and $N_q \to \infty$ for all $q=1,\ldots, Q$. Consequently, the limiting values of some $N_q/N$'s must be 0. The point estimates and variance estimators in \eqref{eqn:estimates} are still well-defined in this regime. We might expect that the asymptotic results in \eqref{eq::basic-asymptotic-inference} still hold because of large $N_q$'s. However, previous theoretical results do not cover this seemingly easy case due to the possibly diverging $H$.

    \item\label{regime:R3} Large $Q$ and small $N_q$'s. In this regime, the number of arms is large but the sample size with each arm is small. Asymptotically, as $N \to \infty$, we have $Q\to \infty$ and $2\le N_q\le \overline{n} $ for some fixed $\overline{n} \ge 2$. This regime is well suited for many factorial experiments (see Example \ref{eg::factorial-design} below), in which the total number of factor combinations can be much larger than the number of replications in each combination \citep[e.g.,][]{mukerjee2006modern, wu2011experiments}. Although the point estimate and variance estimator in \eqref{eqn:estimates} are still well-defined, we do not expect a simple CLT based on the joint asymptotic Normality of ${\hat{Y}}$ due to the small $N_q$'s. Nevertheless, $\widehat{\gamma} = F^\top {\hat{Y}}$, as a linear transformation of ${\hat{Y}}$, can still satisfy the CLT for some choice of $F$. This regime is reminiscence of the so-called {\it proportional asymptotics} in regression analysis, and even there, statistical inference is still not satisfactory in general \citep[e.g.,][]{el2013robust, lei2018asymptotics, el2018can}. Technically, we need to analyze       $ F^\top {\hat{Y}}$ with the dimension of ${\hat{Y}}$ proportional to the sample size under the randomization model. This is a gap in the literature.

    \item\label{regime:R4} Large $Q$ and $N_q = 1 $ for all $q = 1, \ldots, Q$. This regime is much harder than \ref{regime:R3} because the variance estimator in \eqref{eqn:estimates} is not even well defined due to the lack of replications within each arm \citep[e.g.,][]{espinosa2016bayesian}. Therefore, we need to answer two fundamental questions. First, does $ F^\top {\hat{Y}}$ still satisfy the CLT for some $F$? Second, how do we estimate the variance of $ F^\top {\hat{Y}}$? These two questions are the basis for large-sample Wald-type inference in this regime. Neither has been covered by existing results.

    \item\label{regime:R5} Mixture of \ref{regime:R1}--\ref{regime:R4}. In the most general case, it is possible that the number of treatment levels diverges and the group sizes within different treatment arms vary a lot. Theoretically, we can partition the treatment levels into different types corresponding to the four regimes above. Understanding \ref{regime:R5} relies on understanding \ref{regime:R1}--\ref{regime:R4}.    Due to the difficulties in \ref{regime:R1}--\ref{regime:R4} mentioned above,  a rigorous analysis of \ref{regime:R5} requires deeper understanding of the randomization model. This is another gap in the literature.

\end{enumerate}

For descriptive convenience, we define \ref{regime:R1}--\ref{regime:R4} as \textit{nearly uniform designs} and \ref{regime:R5} as \textit{non-uniform designs}, respectively, based on the heterogeneity of the group sizes across treatment arms. Definitions \ref{def:uniform-design} and \ref{def:non-uniform-design} below make the intuition more precise.

\begin{definition}[Nearly uniform design]\label{def:uniform-design} There exists a positive integer $N_0 > 0$ and absolute constants $\underline{c} \le \overline{c}$, such that $ N_q = c_q  {N}_0$ with  $\underline{c}\le c_q\le \overline{c}$, for all $q=1,\ldots, Q$.
\end{definition}

\begin{definition}[Non-uniform design]\label{def:non-uniform-design}
Partition the treatment arms as $\{ 1,\ldots, Q \}= {\mathcal{Q}}_{\textsc{s}}\cup{\mathcal{Q}}_{\textsc{l}}$ with detailed descriptions below.

(i) ${\mathcal{Q}}_{\textsc{l}}$ contains the arms with large sample sizes.  There exists a positive integer $N_0$ and absolute constants $\underline{c} \le \overline{c}$, such that  $N_q = c_q N_0$ with  $\underline{c}\le c_q\le \overline{c}$, for all $q\in{\mathcal{Q}}_{\textsc{l}}$.

(ii) ${\mathcal{Q}}_{\textsc{s}}$ contains the arms with small sample sizes. There exists a fixed integer $\overline{n}  $ such that
    $N_q \le \overline{n} $ for all $q\in{\mathcal{Q}}_{\textsc{s}}$. Further partition ${\mathcal{Q}}_{\textsc{s}}$ as ${\mathcal{Q}}_{\textsc{s}} = {\mathcal{Q}}_{\textsc{u}} \cup {\mathcal{Q}}_{\textsc{r}}$ where
    \begin{itemize}
        \item ${\mathcal{Q}}_{\textsc{r}}$ contains the arms with replications, that is,  $ 2\le N_q \le \overline{n}$ for all $q\in{\mathcal{Q}}_{\textsc{r}}$;
                \item ${\mathcal{Q}}_{\textsc{u}}$ contains the arms without  replications, that is, $N_q = 1$ for all $q\in{\mathcal{Q}}_{\textsc{u}}$.
    \end{itemize}
\end{definition}

For simplicity, we will use $|{\mathcal{Q}}_{\star}|$ and $N_\star = \sum_{q\in{\mathcal{Q}}_\star} N_q$ to denote  the number of arms and the sample size in ${\mathcal{Q}}_\star$, respectively, where $\star \in \{ \textsc{s},\textsc{u},\textsc{r},\textsc{l}\}$. As a special case of Definition \ref{def:non-uniform-design}, $|{\mathcal{Q}}_{\textsc{r}}| = |{\mathcal{Q}}_{\textsc{l}}| = 0$ corresponds to \textit{unreplicated designs} in which each treatment level has only one observation.

We will also use the $2^K$ factorial design as a canonical example for many theoretical results throughout. We review the basic setup of the $2^K$ factorial design in Example \ref{eg::factorial-design}  below \citep{dasgupta2015causal, lu2016covariate, zhao2021regression}.

\begin{example}[Factorial design]
\label{eg::factorial-design}
A $2^K$ factorial design has $K$ binary factors which generate $Q=2^K$ possible treatment levels. Index the potential outcomes $Y_i(q)$'s also as $Y_i(z_1, \ldots, z_K)$'s, where $q = 1,\ldots, Q$ and $z_1,\ldots, z_K = 0,1$. The parameter of interest $    \gamma = F^\top \overline{Y}  $ may consist of a subset of the factorial effects.   The contrast matrix $F$ has orthogonal columns and  $\pm Q^{-1}$ entries; see \citet{dasgupta2015causal} for precise definitions of main effects and interactions.
\end{example}

By definition, the factorial design can either be a nearly uniform design or a non-uniform design, depending on the group sizes across treatment levels. Previous asymptotic results only covered factorial designs under \ref{regime:R1} with fixed $K$ and large sample sizes for all treatment levels. This asymptotic regime can be a poor approximation to finite-sample properties of factorial designs with even a moderate $K$ (for example, if $K=10$ then $Q=2^K > 1000$). Based on simulation, \citet[][Appendix D]{zhao2021regression} showed that CLTs are likely to hold even with diverging $K$ and small sample sizes for all treatment levels.  Allowing for a diverging $K$, \citet[][Theorem A1]{li2017general} derived the CLT for a single factorial effect under the sharp null hypothesis of no treatment effects for any units whatsoever, i.e., $Y_i(1) = \cdots = Y_i(Q)$ for all $i=1,\ldots, N$. However, deriving general asymptotic results for the factorial design has been an open problem in the literature.

\subsection{Our contributions}

Section \ref{sec::open-questions} has reviewed various designs and the associated open problems. In this paper, we will give a unified study of all the designs above. We further the literature in the following ways.

First, we formulate the inference problem under the randomization model in terms of linear permutational statistics. This formulation allows us to build upon the existing results in probability theory \citep{Bolthausen1984AnEO, chatterjee2008multivariate} to derive BEBs on the point estimator of the causal effect. In particular, our analysis emphasizes the dependence on the number of treatment levels and the dimension of the causal effects of interest. Importantly, we derive BEBs that can deal with non-uniform designs with varying group sizes.

Second, we establish a novel BEB on quadratic forms of the linear estimator under the randomization model. Importantly, this BEB allows that the number of treatment levels diverges, the sample sizes across treatment levels vary, and the dimension of the causal effects of interest diverges. It serves as the basis for the $\chi^2$ approximation for large-sample Wald-type inference.

Third, we propose variance estimators for unreplicated designs and mixture designs that allow for the group size to be one in many treatment levels. To the best of our knowledge, the variance estimators are new in the literature of design-based causal inference, although they share some features with those in finely stratified survey sampling \citep[e.g.,][]{cochran1977sampling, wolter2007introduction, breidt2016nonparametric} and experiments \citep[e.g.,][]{abadie2008estimation, fogarty2018mitigating}. However, the theoretical analysis of the new variance estimators is much more challenging because of the dependence of the treatment indicators under the randomization model. We also study their probability limits and establish a complete theory that allows for large-sample Wald-type inference.

Fourth, in the process of achieving the above three sets of results, we established some immediate theoretical results that are potentially useful for other problems. For instance, we prove a novel BEB for linear permutational statistic over convex sets, building upon a recent result based on Stein's method \citep{fang2015rates}.  We also obtain fine results on the sample moments under the randomization model. Due to the space limit, we relegate them to Appendices A and C in the supplementary material.

\subsection{Notation}
We use $C$ to denote generic constants that may vary. Let $\Phi(t)$ denote the cumulative distribution function of a standard Normal distribution. For two sequences of numbers, $a_N$ and $b_N$, let $a_N = O(b_N)$ denote $ a_N \le Cb_N$ for some positive constant $C>0$, and let $a_N = o(b_N) $ denote $a_N/b_N \to 0 $ as $N\to\infty$.
For a  positive integer $N$, define $[N] = \{1,\cdots,N\}$. Let ${\boldsymbol{0}}_N$ and ${\boldsymbol{1}}_N$ denote, respectively, vectors of all zeros and ones in ${\mathbb{R}}^N$. For two random variables $X$ and $X'$, we use $X\lesssim X'$ or $X'\gtrsim X$ to represent that $X'$  stochastically dominates $X$, i.e., ${\mathbb{P}}\{X' \le t\} \le {\mathbb{P}}\{X \le t\}$ for all $t\in{\mathbb{R}}$. For any covariance matrix $V$, let $V^\star$ denote the corresponding correlation matrix.

Consider a matrix $M=(M(h,l))\in{\mathbb{R}}^{H\times H}$. Let $M(\cdot,l)\in {\mathbb{R}}^{H\times 1}$ and $M(h, \cdot)\in{\mathbb{R}}^{1\times H}$ denote its $l$-th column and $h$-th row, respectively.
Let $\varrho_k(M)$ denote its $k$-th largest singular value. Specially, let $\varrho_{\max}(M)$ and $\varrho_{\min}(M)$ denote the largest and smallest singular values, respectively. Define its condition number as the ratio of its largest and smallest singular values: $
\kappa(M) = {\varrho_{\max}(M)}/{\varrho_{\min}(M)}
$.   Let $\|M\|_\textsc{f} = (\sum_{h=1}^H\sum_{l=1}^H M_{hl}^2)^{1/2}$, $\|M\|_{\operatorname{op}} = \{\varrho_{\max}(M^\top M)\}^{1/2}$, $\|M\|_{p,q}= (\sum_{l=1}^H \|M(\cdot,l)\|_p^q)^{1/q} = \{\sum_{l=1}^H(\sum_{h=1}^H   |M_{hl}|^p)^{q/p}\}^{1/q}$ ($1\le p < \infty$), $\|M\|_\infty = \max_{h,l\in [H]} |M_{hl}|$ be, respectively, the Frobenius norm, the operator norm, the $L_{p,q}$ norm and the vectorized $\ell_\infty$ norm.

Design-based results rely crucially on conditions on
$$
M_N(q) = \max_{i\in[N]}|Y_i(q)-\overline{Y}(q)|, \quad (q=1,\ldots, Q)
$$
which is the maximum absolute deviation from the mean for potential outcome $Y_i(q)$'s. \citet{hajek1960limiting} used it in proving the CLT for simple random sampling, and \citet{li2017general} used it in proving CLTs for design-based causal inference. It will also appear frequently in our presentation.

\section{BEBs for the moment estimator under completely randomized experiments}\label{sec:PCLT-projection}
This section presents the BEBs for the moment estimator ${\hat{\gamma}}$ in \eqref{eqn:estimates} under completely randomized experiments.
Section \ref{sec::BEBs-linear-estimators} presents general BEBs for linear projections of ${\hat{\gamma}}$. Sections \ref{sec::BEB-nearly-uniform-design} and \ref{sec::BEB-non-uniform-design} then apply them to derive useful BEBs for nearly uniform and non-uniform designs, respectively.

\subsection{BEBs on the moment estimator}\label{sec::BEBs-linear-estimators}

To simplify the presentation, standardize $ {\hat{\gamma}} = F^\top {\hat{Y}}$:
\begin{align}\label{eqn:tilde-gamma}
    \tilde{\gamma} = V_{\hat{\gamma}}^{-1/2}({\hat{\gamma}} - \gamma) \quad \text{ with } \quad \E{\tilde{\gamma}} = 0 \text{ and } \Var{\tilde{\gamma}} = I_H.
\end{align}
The standardization \eqref{eqn:tilde-gamma} assumes that the covariance matrix $V_{\hat{\gamma}}$ is not singular. We assume it for convenience without loss of generality. When the contrast matrix $F$ has linearly dependent columns and $V_{\hat{\gamma}}$ becomes degenerate, we can focus on a subset of linearly independent columns of $F$.

Our key results are BEBs on linear projections of $    \tilde{\gamma}$. Condition \ref{cond:well-conditioned} below will facilitate the discussion.

\begin{condition}[Non-degenerate covariance matrix of ${\hat{\gamma}}$]\label{cond:well-conditioned}
There exists $\sigma_F \ge 1$ such that
\begin{align}\label{eqn:well-conditioned}
     F^\top  \diag{N_q^{-1}S(q,q)} F \preceq \sigma^2_F F^\top V_{\hat{Y}} F .
\end{align}
\end{condition}

Theorem \ref{thm:be-proj-standard} below gives a general BEB for linear projections of $\tilde{\gamma}$.

\begin{theorem}[BEBs for linear projections of $\tilde{\gamma}$]\label{thm:be-proj-standard}
Assume complete randomization.
(i) There exists a universal constant $C>0$, such that for any $b\in{\mathbb{R}}^H$
with $ \|b\|_2 = 1$, we have
\begin{align*}
\sup_{t\in{\mathbb{R}}}\left|{\mathbb{P}}\{b^\top \tilde{\gamma} \le t\} - \Phi(t)\right| \le C\left\|b^\top V_{\hat{\gamma}}^{-1/2} F^\top\right\|_\infty \cdot
\max_{q\in [Q]} N_q^{-1} M_N(q) .
\end{align*}
(ii) Further assume Condition \ref{cond:well-conditioned}. There exists a universal constant $C>0$, such that
\begin{align}\label{eqn:uniform-be}
    \sup_{b\in{\mathbb{R}}^H,\|b\|_2 =1}\sup_{t\in{\mathbb{R}}}\left|{\mathbb{P}}\{b^\top \tilde{\gamma}  \le t\} - \Phi(t)\right|
    \le C \max_{i\in[N],q\in[Q]} \min\left\{\text{\em I}, \text{\em II}\right\}.
\end{align}
where
\begin{align}\label{eqn::terms1and2-be-proj-standard}
    \text{\em I} = {\sigma_F} \left|\frac{ Y_i(q)-\overline{Y}(q)}{\sqrt{N_q S(q,q)}}\right|,\quad \text{\em II} = \frac{  \|F(q, \cdot)\|_2 \cdot N_q^{-1}|Y_i(q)-\overline{Y}(q)|}{\sqrt{\varrho_{\min}\{F^\top V_{\hat{Y}} F\}}}.
\end{align}
\end{theorem}

The upper bound in Theorem \ref{thm:be-proj-standard}(i) depends on the choice of $b$, whereas the upper bound in Theorem \ref{thm:be-proj-standard}(ii) is uniform over all $b$. The first step of the proof of Theorem \ref{thm:be-proj-standard} is to formulate $\tilde{\gamma}$ as a linear permutational statistic so that we can apply an existing BEB by \citet{Bolthausen1984AnEO}. The second step of the proof is to derive upper bounds in terms of I and II in \eqref{eqn::terms1and2-be-proj-standard} from two different perspectives, which is non-trivial to the best of our knowledge. See Appendix D for technical details. The second step is the key, which makes Theorem \ref{thm:be-proj-standard} applicable in a wide range of designs.

Before discussing the useful consequences of Theorem \ref{thm:be-proj-standard}(ii), we first comment on Condition \ref{cond:well-conditioned}, which plays a key role in deriving the upper bound in Theorem \ref{thm:be-proj-standard}(ii).
Condition \ref{cond:well-conditioned}, coupled with \eqref{eqn:mean-var-Gamma}, implies that
\begin{align*}
    \sigma_F^{-2} F^\top \diag{N_q^{-1}S(q,q)} F \preceq V_{\hat{\gamma}} \preceq F^\top \diag{N_q^{-1}S(q,q)} F,
\end{align*}
i.e., the covariance matrix $V_{\hat{\gamma}}$ is upper and lower bounded by $F^\top \diag{N_q^{-1}S(q,q)} F$, up to constants. Lemma \ref{lem:suff-conds} below gives two sufficient conditions for Condition \ref{cond:well-conditioned} to aid the interpretation.

\begin{lemma}[Sufficient conditions for Condition \ref{cond:well-conditioned}]\label{lem:suff-conds}
(i) Condition \ref{cond:well-conditioned} holds with $\sigma_F = 1$ if the individual causal effects are constant, that is, $\breve{\tau}_i' = F^\top\{(Y_i(q))_{q=1}^Q - \overline{Y}\}= 0$ for all $i\in[N]$.
(ii) Condition \ref{cond:well-conditioned} holds with $\sigma_F = \sigma$ if the condition number of the correlation matrix corresponding to $V_{\hat{Y}}$ is upper bounded by $\sigma^2$.
\end{lemma}

The sufficient conditions in Lemma \ref{lem:suff-conds} are somewhat standard in the literature, especially under \ref{regime:R1}. \citet[][Corollary 2]{li2017general} gives a CLT under the assumption of constant  individual causal effects, which is a special case of Lemma \ref{lem:suff-conds}(i). \citet[][Theorem 5]{li2017general} proves a CLT under the assumption that $S$ has a finite limiting value. When the limit is positive definite, $V_{\hat{Y}}$ also converges to a positive definite matrix, which becomes a special case of  Lemma \ref{lem:suff-conds}(ii).   The general forms of conditions in Condition \ref{cond:well-conditioned} and  Lemma \ref{lem:suff-conds} are more useful for \ref{regime:R2}--\ref{regime:R5}.

Now we comment on Theorem \ref{thm:be-proj-standard}(ii). The upper bound in \eqref{eqn:uniform-be} depends on the minimum value of two terms. It is convenient to apply these two terms to different treatment arms based on the structure of the design. We elaborate this idea by revisiting \ref{regime:R1} to \ref{regime:R5}.
\begin{itemize}
    \item For \ref{regime:R1} and \ref{regime:R2}, because the $N_q$'s are large, we can use term I in \eqref{eqn::terms1and2-be-proj-standard} and obtain a sufficient condition for a vanishing upper bound.

    \item For \ref{regime:R3} and \ref{regime:R4}, the $N_q$'s are bounded and term I in \eqref{eqn::terms1and2-be-proj-standard} has constant order. However, term II in \eqref{eqn::terms1and2-be-proj-standard} is small under mild conditions on $F$. For instance, in the factorial design in Example \ref{eg::factorial-design}, the following algebraic facts hold:
\begin{equation}
\label{eq::contrast-matrix}
\|F\|_{\infty}  = Q^{-1},\quad \|  F(q, \cdot) \|_2 = Q^{-1} \sqrt{H}, \quad \varrho_{\min}(F^\top  F)  = Q^{-1}.
\end{equation}
Combining Condition \ref{cond:well-conditioned} and  \eqref{eq::contrast-matrix}, we have
\begin{align}
    \varrho_{\min}(F^\top V_{\hat{Y}} F) &\ge \sigma_F^{-2} \varrho_{\min}(F^\top \diag{N_q^{-1}S(q,q)} F) \notag\\
    &\ge \sigma_F^{-2}\min_{q\in[Q]} \{N_q^{-1}S(q,q)\}\cdot \varrho_{\min}(F^\top  F) \notag\\
    &= Q^{-1} \sigma_F^{-2} \min_{q\in[Q]} \{N_q^{-1}S(q,q)\}.\label{eqn:factorial-F}
\end{align}
Because the $N_q$'s are bounded, term II in \eqref{eqn::terms1and2-be-proj-standard} has order $O(\sqrt{H/Q})$, which is small if $H/ Q \rightarrow 0$.

\item For \ref{regime:R5}, we can partition the treatment arms based on the sizes of $N_q$'s to achieve a trade-off between terms I and II in \eqref{eqn::terms1and2-be-proj-standard}.  In particular, for non-uniform designs in Definition \ref{def:non-uniform-design}, a natural partition is $ [Q] = {\mathcal{Q}}_{\textsc{l}} \cup {\mathcal{Q}}_{\textsc{S}}$. On the one hand, the arms in ${\mathcal{Q}}_{\textsc{l}}$ contain many units, so term I in \eqref{eqn::terms1and2-be-proj-standard} vanishes asymptotically. On the other hand, ${\mathcal{Q}}_{\textsc{s}}$ contains many arms which makes term II in \eqref{eqn::terms1and2-be-proj-standard} small under mild conditions on $F$.
\end{itemize}

We will provide rigorous results in the next two sections by applying Theorem \ref{thm:be-proj-standard} to obtain useful BEBs for different designs.

\subsection{A BEB with a proper contrast in nearly uniform designs}\label{sec::BEB-nearly-uniform-design}

In \ref{regime:R1} with a fixed $Q$ and large $N_q$'s, it is intuitive to have CLTs for linear transformations of ${\hat{Y}}$. In other regimes, for instance \ref{regime:R4}, CLTs for linear transformations of ${\hat{Y}}$ are less intuitive. Consider a diverging $Q$ and bounded $N_q$'s. If   $F = (1,0,\dots,0)^\top \in {\mathbb{R}}^Q$, then the CLT for $F^\top {\hat{Y}} = {\hat{Y}}_1$ does not hold due to the bounded sample size in treatment arm 1. As another toy example, if $F = ({\boldsymbol{1}}_Q, {\boldsymbol{1}}_Q) \in{\mathbb{R}}^{Q\times 2}$, then $F^\top{\hat{Y}}$ has degenerate covariance structure and Theorem \ref{thm:be-proj-standard} cannot be directly applied. Therefore, CLTs should be established for proper contrast matrices. Corollary \ref{cor:uniform-design-be} below gives a positive result on the BEBs for proper contrasts. We first introduce Condition \ref{condition::proper} below on $F$.

\begin{condition}[Proper contrast]
\label{condition::proper}
The contrast matrix $F$ satisfies $ \|F\|_\infty \le cQ^{-1}$ and $\varrho_{\min}\{F^\top F\} \ge c'Q^{-1}$ for some constants $c,c'>0$.
\end{condition}

Condition \ref{condition::proper} appears to be dependent on the scale of $F$ although the BEB should not depend on the scale of $F$ due tot the standardization of ${\hat{\gamma}}$. We present the above form of Condition \ref{condition::proper} to facilitate the discussion of the factorial design in Example \ref{eg::factorial-design}, in which the scale of $F$ is motivated by scientific questions of interest.
When $Q$ is fixed, Condition \ref{condition::proper} holds if $F$ has full column rank. So in \ref{regime:R1}, Condition \ref{condition::proper} does not impose any additional assumptions beyond the standard ones. When $Q$ diverges, Condition \ref{condition::proper}  rules out sparse $F$ that only results in linear combination of ${\hat{Y}}$ over a small number of treatment arms. Also, the minimum eigenvalue condition in Condition \ref{condition::proper} ensures non-degenerate covariance structure.  Example \ref{eg::factorial-proper-F} below gives more detailed discussion in the factorial design.

 We then give Corollary \ref{cor:uniform-design-be} below.

\begin{corollary}[BEB for nearly uniform designs]\label{cor:uniform-design-be}
Assume complete randomization that satisfies Definition \ref{def:uniform-design}, Conditions \ref{cond:well-conditioned} and \ref{condition::proper}.
There exists a universal constant $C>0$, such that
\begin{align}\label{eqn:uniform-design-be}
 \sup_{b\in{\mathbb{R}}^H,\|b\|_2 =1}\sup_{t\in{\mathbb{R}}}\left|{\mathbb{P}}\{b^\top \tilde{\gamma}  \le t\} - \Phi(t)\right|
 \le
 C\sigma_F    \frac{  \max_{ q\in[Q]} M_N(q)  }{  \{  \min_{q\in[Q]} S(q,q) \}^{1/2} }
\sqrt{ \frac{H}{N} }
 .
\end{align}
\end{corollary}

We make several comments on Corollary \ref{cor:uniform-design-be}. First, technically, we derive the upper bound in Corollary \ref{cor:uniform-design-be} based on the upper bound from term II in \eqref{eqn::terms1and2-be-proj-standard} in Theorem \ref{thm:be-proj-standard}(ii), which is uniform over $b$. Therefore, the upper bound in Corollary \ref{cor:uniform-design-be} preserves the uniformity and does not depend on $b$. Second, the upper bound in \eqref{eqn:uniform-design-be} reveals the interplay of several parameters: the number of contrasts $H$, the number of units $N$, the scale of the potential outcomes $M_N(q)$, the minimum second moments $\min_{q\in[Q]} S(q,q)$  as well as the structure of $F$. Third, the upper bound in \eqref{eqn:uniform-design-be} decreases at the rate of $(H / N)^{1/2} $ which can deal with {\ref{regime:R1}--\ref{regime:R4}}.

\begin{remark}
The denominator of \eqref{eqn:uniform-design-be} depends on $\min_{q\in[Q]} S(q,q)$, which is useful when the variances of the potential outcomes are lower bounded.  For the ease of presentation, we did not discuss more complicated cases such as some $S(q,q)$'s are small. We can slightly modify the proof of Corollary \ref{cor:uniform-design-be} to cover scenarios where some $S(q,q)$'s are close or equal to zero.
\end{remark}

We conclude this section with an example on the nearly uniform factorial design.

\begin{example}[Nearly uniform factorial design]\label{eg::factorial-proper-F}
Recall Example \ref{eg::factorial-design} and assume it satisfies Definition \ref{def:uniform-design}. Let $F\in{\mathbb{R}}^{Q\times H}$ with $H = K+K(K-1)/2={K(K+1)}/{2}$ be the contrast matrix for all main effects and two-way interactions.
Assume Condition \ref{cond:well-conditioned} and recall \eqref{eq::contrast-matrix}.  Corollary \ref{cor:uniform-design-be} implies
\begin{align} \label{eqn:factorial-BE}
 \sup_{b\in{\mathbb{R}}^H,\|b\|_2 =1}\sup_{t\in{\mathbb{R}}}\left|{\mathbb{P}}\{b^\top \tilde{\gamma}  \le t\} - \Phi(t)\right|
 \le
 C\sigma_F     \frac{ \max_{q\in[Q]} M_N(q) }{ \{  \min_{q\in[Q]} S(q,q) \}^{1/2} }  \sqrt{\frac{K^2}{N}}   .
\end{align}
From \eqref{eqn:factorial-BE}, we can obtain a sufficient condition for the the upper bound to converge to 0, which implies a CLT of $\tilde{\gamma} $.
\end{example}

\subsection{A BEB for non-uniform designs}\label{sec::BEB-non-uniform-design}

Now consider the non-uniform design in Definition \ref{def:non-uniform-design}.
We can apply Theorem \ref{thm:be-proj-standard} to establish the following BEB for non-uniform designs:

\begin{corollary}[BEB for non-uniform designs]\label{cor:non-uniform-design-be}
Assume complete randomization  that satisfies Definition \ref{def:non-uniform-design} and Conditions \ref{cond:well-conditioned} and \ref{condition::proper}. Also assume
\begin{align}\label{eqn:small-H}
    Q \ge 2(c^2/c')H|{\mathcal{Q}}_{\textsc{l}}|
\end{align}
recalling the constants $c$ and $c'$ in Condition \ref{condition::proper}.
There exists a universal constant $C>0$, such that
\begin{align}\label{eqn:non-uniform-design-be}
 &\sup_{b\in{\mathbb{R}}^H,\|b\|_2 =1}\sup_{t\in{\mathbb{R}}}\left|{\mathbb{P}}\{b^\top \tilde{\gamma}  \le t\} - \Phi(t)\right|  \\
 \le&  C \sigma_F  \max\left\{ \max_{q\in{\mathcal{Q}}_{\textsc{l}}}\frac{  M_N(q)}{\sqrt{N_q S(q,q)}}, \frac{   \max_{q\in {\mathcal{Q}}_{\textsc{s}}} M_N(q)}{   \{   \min_{q\in {\mathcal{Q}}_{\textsc{s}}} S(q,q) \}^{1/2}  }\cdot \sqrt{\frac{H}{N_{\textsc{s}}}}\right\}. \notag
\end{align}
\end{corollary}

There are two key steps in applying Theorem \ref{thm:be-proj-standard} to prove Corollary \ref{cor:non-uniform-design-be}. The first step is to partition ${\mathcal{Q}}$ into ${\mathcal{Q}}_{\textsc{s}} \cup {\mathcal{Q}}_{\textsc{l}}$  based on the size of the arms.  The second step is to simplify the denominator of term II of \eqref{eqn:uniform-be}. The obtained upper bound  \eqref{eqn:non-uniform-design-be} is uniform over all $b$. Moreover, it depends on the sizes of the treatment arms in a subtle way.  On the one hand, for $q\in{\mathcal{Q}}_{\textsc{l}}$, the $N_q$'s are large, so the first part of \eqref{eqn:non-uniform-design-be} converges to zero if the following ``local'' condition holds for all $q\in{\mathcal{Q}}_{\textsc{l}}$:
$$
\frac{  M_N(q)^2  }{ S(q,q) } = o(N_q)  .
$$
On the other hand, for $q\in{\mathcal{Q}}_{\textsc{s}}$, the $N_q$'s are small, but the second part of \eqref{eqn:non-uniform-design-be} still converges to zero if the following ``global'' condition holds:
$$
\frac{   \max_{q\in {\mathcal{Q}}_{\textsc{s}}} M_N(q)^2 }{     \min_{q\in {\mathcal{Q}}_{\textsc{s}}} S(q,q)   } = o\left(  \frac{H}{ N_{\textsc{s}} } \right).
$$

We conclude this section with an example on the non-uniform factorial design.

\begin{example}[Non-uniform factorial design]
Recall Example \ref{eg::factorial-design}.  Assume the baseline arm $q=1$ contains a large number of units possibly due to lower cost while the other arms have $N_q \le \overline{n}$ for some fixed $\overline{n}$. This gives a non-uniform design by Definition \ref{def:non-uniform-design} with ${\mathcal{Q}}_{\textsc{l}} = \{1\}$ and ${\mathcal{Q}}_{\textsc{s}} = \{2, \ldots, Q\}$.
Let $F\in{\mathbb{R}}^{Q\times H}$ with $H={K(K+1)}/{2}$ be the contrast matrix for all main effects as well as two-way interactions.
Assume Condition \ref{cond:well-conditioned} and recall \eqref{eq::contrast-matrix}. Applying Corollary \ref{cor:non-uniform-design-be}, we have
\begin{align}\label{eqn:non-uniform-factorial-be}
 &\sup_{b\in{\mathbb{R}}^H,\|b\|_2 =1}\sup_{t\in{\mathbb{R}}}\left|{\mathbb{P}}\{b^\top \tilde{\gamma}  \le t\} - \Phi(t)\right|  \\
 \le&  C \sigma_F  \max\left\{\frac{ M_N(1)}{\sqrt{N_1 S(1,1)}},
 \frac{  \max_{ q \geq 2}   M_N(q)}{    \{  \min_{q\geq 2} S(q,q) \}^{1/2} }  \sqrt{  \frac{K^2}{N_\textsc{s}} }  \right\}. \notag
\end{align}
From \eqref{eqn:non-uniform-factorial-be}, if $K\to\infty, N_1\to\infty$, and
$$
\frac{M_N(1)}{\sqrt{S(1,1)}} = o(N_1^{1/2}) ,\quad
\frac{  \max_{q\geq 2} M_N(q)}{ \{  \min_{q\geq 2} S(q,q)\}^{1/2} } = o( {  N_\textsc{s}^{1/2} } / K),
$$
then the upper bound in \eqref{eqn:non-uniform-factorial-be} vanishes asymptotically.
\end{example}

\section{Design-based causal inference}\label{sec:inference}

Now we turn to the central task of design-based causal inference under complete randomization. We focus on the large-sample Wald-type inference based on the quadratic form
$$
\widehat{T} = ({\hat{\gamma}} - \gamma)^\top {\hat{V}}_{{\hat{\gamma}}}^{-1} ({\hat{\gamma}} - \gamma),
$$
recalling the point estimator  ${\hat{\gamma}}$ and  the variance estimator ${\hat{V}}_{{\hat{\gamma}}}$ in \eqref{eqn:estimates}. In \ref{regime:R1} with fixed $(Q, H)$ and large $N_q$'s, the standard asymptotic argument suggests that we can use $q_\alpha$, the upper $\alpha$-quantile of $\chi^2(H)$, as the critical value for the  the quadratic form. For simplicity, we say that the corresponding confidence set is asymptotic valid if  $ \lim_{N\to\infty} \mathbb{P}\{  \widehat{T}  \le q_\alpha \} \ge 1-\alpha$.

The rigorous theoretical justification for the above Wald-type inference procedure typically follows from two steps:
\begin{enumerate}
\item
[(Step 1)]\label{step1}
first analyze the asymptotic distribution of the corresponding quadratic form with the true covariance matrix
\begin{align}\label{eqn:quad-form}
   T = ({\hat{\gamma}} - \gamma)^\top V_{{\hat{\gamma}}}^{-1} ({\hat{\gamma}} - \gamma).
\end{align}

\item
[(Step 2)]\label{step2}
second construct a consistent or conservative estimator ${\hat{V}}_{{\hat{\gamma}}}$ for the true covariance matrix    $V_{{\hat{\gamma}}}$.
\end{enumerate}

Under regime \ref{regime:R1}, both Steps 1 and 2 have rigorous theoretical justification ensured by \eqref{eq::basic-asymptotic-inference}. Beyond \ref{regime:R1}, it is challenging to derive the asymptotic distribution of the quadratic form in Step 1 especially when $H$ and thus the degrees of freedom of $T$ diverge. To achieve the requirement in Step 1, we use results based on Stein's method to derive BEBs on quadratic forms of linear permutational statistics. To avoid excessive notation, we present the results that are most relevant to our inference problem in the main paper and relegate more general yet more complicated results to Appendices A and C. Moreover, the sample variances ${\hat{S}}(q,q)$'s and thus the variance estimator ${\hat{V}}_{{\hat{\gamma}}}$ in \eqref{eqn:estimates} are not even well defined when some treatment arms do not have replications of the outcome. Without replications in all arms, we must find an alternative form of  ${\hat{V}}_{{\hat{\gamma}}}$ to estimate $V_{{\hat{\gamma}}}$. This is a salient problem for \ref{regime:R4} and \ref{regime:R5}. Finally, in all regimes \ref{regime:R1}--\ref{regime:R5}, we need to study the properties of ${\hat{V}}_{{\hat{\gamma}}}$ to achieve the requirement in Step 2.

Due to the different levels of technical complexities, we divide this section into three subsections.   Section \ref{sec:var-uniform} discusses nearly uniform designs with replications in all arms. Section \ref{sec:var-unreplicate} discusses unreplicated designs. Section \ref{sec::non-uniform-design-inference} discusses the general non-uniform designs. In every subsection, we first present a BEB on the quadratic form in Step 1, then present the properties of the covariance estimator ${\hat{V}}_{{\hat{\gamma}}}$, and finally present the formal result to justify Wald-type inference.

To facilitate the discussion, we introduce the following notation
$$
T_0 = \xi_H^\top \xi_H \quad \text{ where }
\xi_H\sim{\mathcal{N}}(0,I_H)
$$
for a $\chi^2(H)$ random variable with possibly diverging dimension. It has mean $H$ and variance $2H$. Both $\widehat{T}$ and $T$ are related to $T_0.$ Moreover, we introduce the following moment condition on the potential outcomes.

 \begin{condition}[Bounded fourth moment of the potential outcomes]\label{cond:moments}
There exists a $\Delta>0$ such that
$
    \max_{q\in[Q]}  N^{-1} \sum_{i=1}^N \{Y_i(q) - \overline{Y}(q)\}^4 \le \Delta^4.
$
\end{condition}

\subsection{Nearly uniform design with replications in all arms}\label{sec:var-uniform}

In this subsection, we study the Wald-type inference for nearly uniform designs given by Definition \ref{def:uniform-design}.
First, we present a BEB for $T$ in \eqref{eqn:quad-form} in Theorem \ref{thm:quad-be-nearly-uniform} below.

\begin{theorem}[BEB for  the quadratic form $T$ for  nearly uniform designs with replications]\label{thm:quad-be-nearly-uniform}
Assume complete randomization that satisfies Definition \ref{def:uniform-design}. Assume Conditions \ref{cond:well-conditioned} and \ref{condition::proper}. There exists a universal constant $C>0$, such that
\begin{align}\label{eqn:quad-be-nearly-uniform}
    \sup_{t\in{\mathbb{R}}} |{\mathbb{P}}(T\le t)-{\mathbb{P}}(T_0 \le t)| \le \frac{C\max_{q\in[Q]}M_N(q)^3 }{ \{ \min_{q\in[Q]}S(q,q) \}^{3/2}}\cdot \frac{H^{19/4}}{N^{1/2}}.
\end{align}
\end{theorem}

Theorem \ref{thm:quad-be-nearly-uniform} bounds the difference between $T$ and $T_0$ with possibly diverging $H$. Its upper bound is more useful when $H^{19/2} / N \rightarrow 0$, which restricts the number of contrasts of interest. The condition $H^{19/2} / N \rightarrow 0$ holds naturally in the factorial design in Example \ref{eg::factorial-design} under regime \ref{regime:R4} if only the main effects and two-way interactions are of interest with $H = O(K^2) = O( (\log N)^2 ) $.

Second, we discuss variance estimation. Recall ${\hat{S}}(q,q)$ and ${\hat{V}}_{\hat{Y}}$ be defined as in  \eqref{eqn:hS} and \eqref{eqn:hY-hV}. Consider the point estimator ${\hat{\gamma}}$ and covariance estimator ${\hat{V}}_{\hat{\gamma}}$ in \eqref{eqn:estimates}.
We have  Theorem \ref{thm:uniform-var} below.

\begin{theorem}[Variance estimation in nearly uniform designs]\label{thm:uniform-var}
Consider designs that satisfies Definition \ref{def:uniform-design} with $\min_{q\in[Q]} N_q \ge 2$. Assume Condition \ref{cond:moments}.
\begin{enumerate}[label = (\roman*)]
    \item\label{thm:uniform-var-1}
${\mathbb{E}}\{{\hat{V}}_{{\hat{\gamma}}}\} \succeq  V_{\hat{\gamma}}$.

 \item\label{thm:uniform-var-2}
 $ \|{\hat{V}}_{{\hat{\gamma}}}-{\mathbb{E}}\{{\hat{V}}_{{\hat{\gamma}}}\}\|^2_{\infty} = O_{{\mathbb{P}}}\left({\|F\|_\infty^4 Q^4 N^{-3}\Delta^4 H^2}\right).$

 \item \label{thm:uniform-var-3}
 $ \|{\hat{V}}_{{\hat{\gamma}}}-{\mathbb{E}}\{{\hat{V}}_{{\hat{\gamma}}}\}\|^2_{{\operatorname{op}}} = O_{{\mathbb{P}}}\left({ \|F\|_\infty^4 Q^4 N^{-3}\Delta^4 H^4}\right).$

\end{enumerate}

\end{theorem}

Theorem \ref{thm:uniform-var}\ref{thm:uniform-var-1} reviews that the covariance estimator ${\hat{V}}_{{\hat{\gamma}}}$ is conservative, which is well-known in design-based causal inference \citep{neyman1923application, imbens15, li2017general}. Theorem \ref{thm:uniform-var}\ref{thm:uniform-var-2} and \ref{thm:uniform-var-3} are novel results on the stochastic orders of of the estimation error of the covariance estimator in $L_{\infty}$ norm and operator norm, respectively. In Example \ref{eg::factorial-design} of the factorial design with $\|F\|_\infty = O(Q^{-1})$, if $\Delta$ is constant, then Theorem \ref{thm:uniform-var} simplifies to
$$
   N\|{\hat{V}}_{{\hat{\gamma}}}-{\mathbb{E}}\{{\hat{V}}_{{\hat{\gamma}}}\}\|_{\infty} =O_{{\mathbb{P}}}\left(  H/N^{1/2}  \right),\quad
    N\|{\hat{V}}_{{\hat{\gamma}}}-{\mathbb{E}}\{{\hat{V}}_{{\hat{\gamma}}}\}\|_{\operatorname{op}}  = O_{{\mathbb{P}}}\left( H^2/N^{1/2}   \right).
$$
The estimation error shrinks to zero quickly if only the main effects and two-way interactions are of interest.
The results in Theorem \ref{thm:uniform-var} suffice for inference, and we relegate the finer probability tail bound for  ${\hat{V}}_{{\hat{\gamma}}}$ to the supplementary material.

Third, we present formal results on inference. To simplify the presentation, we impose Condition \ref{cond:easy-spec} below.

\begin{condition}\label{cond:easy-spec}
(i) The absolute value of the centered potential outcomes are upper bounded by a constant $\nu>0$.
(ii) The variance of the potential outcomes are lower bounded a constant $\underline{S} > 0$.
\end{condition}

The boundedness in potential outcomes in Condition \ref{cond:easy-spec} is not necessary but convenient for the theory. We can relax it to allow for light tails of the potential outcomes. We omit the results due to the technicalities.  Theorem \ref{thm:wald-uniform} below justifies the Wald-type inference under the nearly uniform design with replications, where $H$ can be either fixed or diverging.

\begin{theorem}[Wald-type inference under replicated nearly uniform design]\label{thm:wald-uniform}
Consider the nearly uniform design given by Definition \ref{def:uniform-design} with $\min_{q\in[Q]}N_q\ge 2$. Assume Conditions \ref{cond:well-conditioned}-\ref{cond:easy-spec}.
Define
$
    W_N = V_{\hat{\gamma}}^{1/2}\E{{\hat{V}}_{\hat{\gamma}}}^{-1}V_{\hat{\gamma}}^{1/2} \in {\mathbb{R}}^{H\times H}.
$
\begin{enumerate}[label = (\roman*)]
    \item For a fixed $H$, assume there exists a {$W_\infty\in{\mathbb{R}}^{H\times H}$} such that
    $ \lim_{N\to\infty}W_N = W_\infty $.
    Use ${\mathcal{L}}$ to denote the distribution of
    $
    \xi_H^\top W_\infty \xi_H , \text{ where } \xi_H \sim {\mathcal{N}}_H(0,I_H).
    $
    We have
    \begin{align*}
        ({\hat{\gamma}} - \gamma )^\top {\hat{V}}_{{\hat{\gamma}}}^{-1} ({\hat{\gamma}} - \gamma) \rightsquigarrow {\mathcal{L}}  , ~N\to\infty.
    \end{align*}
    and ${\mathcal{L}} \lesssim \chi^2(H)$.
    Moreover, the Wald-type confidence set is asymptotically valid.

    \item For a diverging $H$,     we have
    \begin{align*}
        \frac{({\hat{\gamma}} - \gamma)^\top  {\hat{V}}_{{\hat{\gamma}}}^{-1} ({\hat{\gamma}} - \gamma) - \trace{W_N}}{\sqrt{2\trace{W_N^2}}} \rightsquigarrow {\mathcal{N}}(0,1) \text{ if } ~H,N\to\infty,~ H^{19/4}N^{-1/2} \to 0.
    \end{align*}
    Moreover, the Wald-type confidence set is asymptotically valid.
\end{enumerate}
\end{theorem}

Theorem \ref{thm:wald-uniform}(i) reviews the known result for \ref{regime:R1} with a fixed $H$. Theorem \ref{thm:wald-uniform}(ii) is a novel result that allows for diverging $Q$ and $H$. It relies crucially on the BEB on the quadratic form $T$ in Theorem \ref{thm:quad-be-nearly-uniform} and the stochastic properties of ${\hat{V}}_{\hat{\gamma}}$ in Theorem \ref{thm:uniform-var}.

\subsection{Unreplicated design}\label{sec:var-unreplicate}

In this subsection, we study inference for unreplicated designs with $N_q = 1$ for $q=1,\ldots, Q$. On the one hand, the BEB on $T$ is identical to that in Theorem \ref{thm:quad-be-nearly-uniform}. We give the formal result in  Theorem \ref{thm:quad-be-unreplicated} below for completeness.

\begin{theorem}[BEB for the quadratic form in unreplicated designs]\label{thm:quad-be-unreplicated}
Assume complete randomization with $N_q = 1$ for $q=1,\ldots, Q$. Assume Conditions \ref{cond:well-conditioned} and \ref{condition::proper}. The BEB \eqref{eqn:quad-be-nearly-uniform} holds.
\end{theorem}

On the other hand, covariance estimation without replications is a fundamentally challenging problem as reviewed in Section \ref{sec:motivation}. The commonly-used  covariance estimator ${\hat{V}}_{\hat{\gamma}}$ in \eqref{eqn:estimates} is not well defined. We must construct a new estimator.  In unreplicated designs, the observed allocation $Z_i$ and the arm $q$ has a one-to-one correspondence. Hence we can denote the single observed outcome in arm $ q$ by $Y_q$. The point estimator still has the form ${\hat{\gamma}} = F^\top {\hat{Y}}$ where ${\hat{Y}} = (Y_1,\ldots, Y_Q)^\top$ is simply the observed outcome vector. Without replications, we cannot calculate ${\hat{S}}(q,q)$ based on only the single observation within arm $q \in {\mathcal{Q}}_\textsc{u} = \{ 1,\ldots, Q\} $. With a little abuse of notation, we still consider the covariance estimator of the form:
\begin{align}\label{eqn:hV-hY-abuse}
    {\hat{V}}_{{\hat{\gamma}}} =  F^\top {\hat{V}}_{\hat{Y}}  F,
\end{align}
where ${\hat{V}}_{\hat{Y}}$ is a $Q\times Q$ diagonal matrix. The key is to construct its diagonal elements ${\hat{V}}_{\hat{Y}}(q,q)$ for all $q$'s.

To obtain substitutes for ${\hat{S}}(q,q)$, we must borrow information across treatment arms. This motivates us to consider the following grouping strategy.

\begin{definition}[Grouping]
\label{def::Grouping}
Partition ${\mathcal{Q}}_\textsc{u}$ as $   {\mathcal{Q}}_{\textsc{u}} = \cup_{g=1}^G {\mathcal{Q}}_{\textsc{u},g}$ where ${\mathcal{Q}}_{\textsc{u},g} \cap {\mathcal{Q}}_{\textsc{u},g'}  = \varnothing$ for all $g\neq g'$ and $ | {\mathcal{Q}}_{\textsc{u},g}| \ge 2$ for all $g \in [G]$. The partition does not depend on the observed data.
\end{definition}

Definition \ref{def::Grouping} does not allow for data-dependent grouping, which can cause theoretical complications. Examples \ref{exp:pairing} and \ref{exp:regression} below are special cases of Definition \ref{def::Grouping}. By the construction in Definition \ref{def::Grouping}, the ${\mathcal{Q}}_{\textsc{u},g}$'s have no overlap, so we can also use $\langle g \rangle$  to denote ${\mathcal{Q}}_{\textsc{u},g}$ and ${\mathcal{G}} = \{ \langle g \rangle\}_{g=1}^G$ to denote the grouping strategy without causing potential confusions.  Moreover, $|\langle g \rangle|$ must be larger than or equal to two so that there are at least two treatment levels in each $\langle g \rangle$. In general, we use $ \langle g \rangle_q$ to indicate the group $\langle g \rangle$ that contains arm $q$, but when no confusions would arise, we also simplify the notation to $ \langle g \rangle$ if the corresponding $q$ is clear from the context.

Define
\begin{align*}
    {\hat{Y}}_{\langle g \rangle} = \frac{1}{|\langle g \rangle|} \sum_{q\in\langle g \rangle} Y_q,
\end{align*}
 as the group-specific average, and construct
 \begin{align}\label{eqn:hV-QU}
{\hat{V}}_{\hat{Y}}(q,q) = \mu_{\langle g \rangle}(Y_q - {\hat{Y}}_{{\langle g \rangle}})^2 ,\quad \text{ if } q\in{\langle g \rangle}
\end{align}
as the $q$th diagonal element of ${\hat{V}}_{\hat{Y}}$, where
\begin{align}\label{eqn:mu-bg}
     \mu_{\langle g \rangle} = (1-2N^{-1})^{-1}(1-|{\langle g \rangle}|^{-1})^{-2}
\end{align}
is a correction factor that is motivated by the theory below. Although the mean of $    {\hat{Y}}_{\langle g \rangle} $ has a simple formula, the mean of ${\hat{V}}_{\hat{Y}}(q,q) $ has a cumbersome form. We present an lower bound of $    {\mathbb{E}}\{{\hat{V}}_{\hat{Y}}(q,q)\} $ below and relegate the complete formula to the supplementary material. The results require Condition \ref{cond:cond-N} below on the largest eigenvalue of the population correlation of the potential outcomes in group $\langle g \rangle$, defined as
\begin{align}\label{eqn:rho-max-S}
    \varrho_{\langle g \rangle} = \varrho_{\max}\{(S^\star(q,q'))_{q,q'\in\langle g\rangle}\} .
\end{align}

\begin{condition}[Bound on $  \varrho_{\langle g \rangle} $]\label{cond:cond-N}
$
    N- \varrho_{\langle g \rangle} -  (|{\langle g \rangle}|-1) \ge 0
$
for all $g\in{\mathcal{G}}$.
\end{condition}

Condition \ref{cond:cond-N} reflects a trade-off between $N$, $\langle g\rangle$ and $\varrho_{\langle g \rangle}$. It is more likely to hold with smaller correlations between arms within the same group and smaller subgroup sizes. Using a natural bound $\rho_{\langle g \rangle} \le |{\langle g \rangle}|$, Condition \ref{cond:cond-N} holds if $|{\langle g \rangle}| \le {(N + 1)}/{2}$ for all $g\in[G]$. Examples \ref{exp:pairing} and \ref{exp:regression} below satisfy Condition \ref{cond:cond-N} automatically.
With Condition \ref{cond:cond-N}, we can present Lemma \ref{lemma::mean-variance-grouping} below.

\begin{lemma}
[Sample mean and variance under grouping]
\label{lemma::mean-variance-grouping}
Assume grouping ${\mathcal{G}}$. We have
\begin{align*}
\E{{\hat{Y}}_{{\langle g \rangle}}}  =    \overline{Y}_{\langle g \rangle} ,\quad \text{ where }   \overline{Y}_{\langle g \rangle} = \frac{1}{|{\langle g \rangle}|N} \sum_{q\in{\langle g \rangle}} \sum_{i=1}^NY_i(q) = \frac{1}{|{\langle g \rangle}|} \sum_{q\in{\langle g \rangle}} \overline{Y}(q).
\end{align*}
Further assume Condition \ref{cond:cond-N}. We have
\begin{align*}
     {\mathbb{E}}\{{\hat{V}}_{\hat{Y}}(q,q)\} \ge S(q,q) + \underbrace{\underline{\Omega}(q,q)}_{\text{\em term III}} +  \underbrace{\mu_{\langle g \rangle}  (\overline{Y}(q) - \overline{Y}_{\langle g \rangle})^2}_{\text{\em term IV}},
\end{align*}
where
\begin{align} \label{eq::Omega-qq}
    {\Omega}(q,q) &= \mu_{\langle g \rangle}|g|^{-2}\left(1-\frac{\varrho_{\langle g \rangle}}{N} - \frac{|g|-1}{N}\right)\sum_{q'\in{\langle g \rangle},q'\neq q}S(q',q')  \geq 0.
\end{align}
\end{lemma}

By Lemma \ref{lemma::mean-variance-grouping}, ${\hat{V}}_{\hat{Y}}(q,q)$, as an estimator for $S(q,q) $, is conservative, and the conservativeness depends on the variation of other arms $q'$ that belong to ${\langle g \rangle}_q$ (term III) and the between-arm heterogeneity in means within ${\langle g \rangle}_q$ (term IV). We comment on some special cases below.
\begin{itemize}
    \item If we assume homogeneity in means within subgroups, i.e.,
\begin{align}\label{eqn:weak-means}
    \overline{Y}(q) = \overline{Y}_{\langle g \rangle},~ \text{ for all } q\in\langle g \rangle,
\end{align}
then term IV vanishes.

   \item If we assume homoskedasticity across treatment arms within the same subgroup, i.e.,
\begin{align}\label{eqn:weak-vars}
    S(q,q) = S(q',q'),~ \text{ for all } q,q'\in \langle g \rangle,
\end{align}
then term III becomes
\begin{align}\label{eqn:strong-po}
    {\Omega}(q,q) = \mu_{\langle g \rangle}(|g|-1)|g|^{-2}\left(1-\frac{\varrho_{\langle g \rangle}}{N} - \frac{|g|-1}{N}\right)S(q,q) .
\end{align}
Then we can combine \eqref{eqn:strong-po} with $S(q,q)$ and use a smaller correction factor
\begin{align*}
    \mu'_{\langle g \rangle} = (1-|g|^{-1})^{-1}\{(1-|g|^{-1})(1-2N^{-1}) + |g|^{-1}(1 - (2|g| - 1)/N) \}^{-1} \le \mu_{\langle g \rangle}
\end{align*}
to reduce the conservativeness of variance estimation.

  \item If we assume the strong null hypothesis within subgroups, i.e.,
  \begin{align*}
      Y_i(q) = Y_i(q'), \text{ for all }i\in[N] \text{ and } q,q'\in {\langle g \rangle},
  \end{align*}
  then both \eqref{eqn:weak-means} and \eqref{eqn:weak-vars} hold.  Applying the correction factor $\mu'_{\langle g \rangle}$, we can show
$
      {\mathbb{E}}\{{\hat{V}}_{\hat{Y}}(q,q)\} = S(q,q).
$
\end{itemize}

Lemma \ref{lemma::mean-variance-grouping} suggests that ideally, we should group treatment arms based on the prior knowledge on the variation and means of the potential outcomes. While more general grouping strategies are possible, we give two examples for their simplicity of implementation. Both target the factorial design in Example \ref{eg::factorial-design}.

\begin{example}[Pairing by the lexicographic order]\label{exp:pairing}
Recall Example \ref{eg::factorial-design}.
We order the observations based on the lexicographical order of their treatment levels, then group the $(2k-1)$-th level with the $(2k)$-th level $(1\le k\le 2^{K-1})$. When $K=3$, the grouping reduces to
 \begin{align*}
   \langle 1 \rangle=  \{ (000), (001) \}, \quad
      \langle 2 \rangle   = \{ (010), (011) \}, \quad
        \langle 3 \rangle   = \{    (100), (101) \} , \quad
          \langle 4 \rangle   =    \{   (110), (111)\} .
 \end{align*}
 If the last factor has small effect on the outcome, then we expect small differences of the mean potential outcomes within grouping.
As a sanity check, Condition \ref{cond:cond-N} holds under this grouping strategy.
\end{example}

\begin{example}[Grouping based on a subset of the factors]\label{exp:regression}
Recall Example \ref{eg::factorial-design} again. If we have the prior knowledge that $K_0 < K$ factors are the most important ones, we can group the treatment levels based on these factors. Without loss of generality, assume that the first $K_0$ factors are the important ones. In particular, we can create $ G = 2^{K_0} < Q$ groups, with each group $     {\langle g \rangle}$ corresponding to treatment levels with the same important factors. Example \ref{exp:pairing} above is a special with the first $K-1$ factors as the important ones.
Also, Condition \ref{cond:cond-N} holds under this grouping strategy.
\end{example}

Now we turn to theoretical analysis of \eqref{eqn:hV-QU}.  Its properties depend on how successful the grouping ${\mathcal{G}}$ is, quantified by Condition \ref{cond:bounded-bgv} below.

\begin{condition}[Bound on the within-group variation in potential outcome means]\label{cond:bounded-bgv}
There exists a $\zeta>0$, such that
$
    \max_{g\in [G]} \max_{q\in{\langle g \rangle}} |\overline{Y}(q) - \overline{Y}_{\langle g \rangle}| \le \zeta.
$
\end{condition}

The $\zeta$ in Condition \ref{cond:bounded-bgv} bounds the between-arm distance of the mean potential outcomes under grouping ${\mathcal{G}}$. It plays a key role in Theorem \ref{thm:unreplicated-var} below.

\begin{theorem}[Variance estimation for unreplicated designs]\label{thm:unreplicated-var}
Consider designs that satisfy Definition \ref{def:non-uniform-design} and the covariance estimator in \eqref{eqn:hV-QU}.
\begin{enumerate}[label = (\roman*)]
    \item\label{thm:unreplicated-var-1}
Assume Condition \ref{cond:cond-N}.    We have

    \begin{align*}
        {\mathbb{E}}\{{\hat{V}}_{\hat{Y}}\}
        =  V_{\hat{Y}}  +  \Omega  +  \diag{\mu_{\langle g \rangle}(\overline{Y}(q) - \overline{Y}_{\langle g \rangle})^2}_{q\in{\mathcal{Q}}_{\textsc{u}}} + N^{-1} (\Theta + S)
    \end{align*}
    with $\Omega = \diag{ \Omega(q,q)  }_{q\in{\mathcal{Q}}_{\textsc{u}}}$ and $ \Theta = \diag{\Theta(q,q)}_{q\in{\mathcal{Q}}_\textsc{U}}$, where the $\Omega(q,q) $'s are defined in \eqref{eq::Omega-qq} and the $\Theta(q,q)$'s are bounded by $0\le \Theta(q,q) \le 5\mu_{\langle g \rangle}\max_{q'\in \langle g \rangle} S(q',q').$
    Therefore,     ${\mathbb{E}}\{F^\top{\hat{V}}_{\hat{Y}} F\} \succeq V_{\hat{\gamma}} $.

    \item \label{thm:var-unreplicated-2}
Assume Conditions \ref{cond:moments} and    \ref{cond:bounded-bgv}. We have

    \begin{align*}
        \|{\hat{V}}_{{\hat{\gamma}}}-{\mathbb{E}}\{{\hat{V}}_{{\hat{\gamma}}}\}\|^2_{\infty} & = O_{\mathbb{P}}\left\{(\max_{g\in[G]}\mu_{\langle g \rangle})^2 \|F\|_\infty^4\Delta^2
    (\Delta^2+\zeta^2)NH^2\right\}.
    \end{align*}
    \item \label{thm:var-unreplicated-3}
    Assume Conditions \ref{cond:moments} and    \ref{cond:bounded-bgv}. We have

    \begin{align*}
        \|{\hat{V}}_{{\hat{\gamma}}}-{\mathbb{E}}\{{\hat{V}}_{{\hat{\gamma}}}\}\|^2_{\operatorname{op}} & = O_{\mathbb{P}}\left\{(\max_{g\in[G]}\mu_{\langle g \rangle})^2 \|F\|_\infty^4\Delta^2
    (\Delta^2+\zeta^2)NH^4\right\}.
    \end{align*}
\end{enumerate}

\end{theorem}

Theorem \ref{thm:unreplicated-var}\ref{thm:unreplicated-var-1} demonstrates that based on $\eqref{eqn:hV-QU}$, the covariance estimator  ${\hat{V}}_{\hat{Y}}$ is conservative for $V_{\hat{Y}}$, which implies that $F^\top{\hat{V}}_{\hat{Y}} F$ is conservative for the true covariance matrix of ${\hat{\gamma}}$. The conservativeness, however, has a totally different pattern compared to the setting with replications within all arms \citep{neyman1923application, imbens15, li2017general}.
Theorem \ref{thm:unreplicated-var}(i) shows three sources of conservativeness. The first part, captured by $\Omega $, is due to the between-arm heteroskedasticity within each subgroup. However, it is fundamentally difficult to estimate each $S(q,q)$ without replications. The second part, captured by $\diag{\mu_{\langle g \rangle}(\overline{Y}(q) - \overline{Y}_{\langle g \rangle})^2}_{q\in[Q]}$,  is due to the between-arm heterogeneity in means within each subgroup.
If the grouping strategy ensures that the grouped arms have similar population average of potential outcomes, the part will be small.  The third part, captured by $N^{-1} (\Theta + S)$,  is due to the difficulty of estimating $S$ and in particular, the off-diagonal terms of $S$. The difficulty of estimating $S$ has been well documented ever since \citet{neyman1923application} even in experiments with replications in each arm. It is possible to reduce this part but it requires additional assumptions, for example, the individual causal effects are constant.

Theorem \ref{thm:unreplicated-var}\ref{thm:var-unreplicated-2} and \ref{thm:var-unreplicated-3} give the stochastic order of the estimation error of the covariance estimator ${\hat{V}}_{{\hat{\gamma}}}$ under the $L_\infty$ norm and operator norm, respectively. If $\max_{g\in[G]}\mu_{\langle g \rangle}$, $\Delta$ and $\zeta$ are all constants, then
$$
    N\|{\hat{V}}_{{\hat{\gamma}}}-{\mathbb{E}}\{{\hat{V}}_{{\hat{\gamma}}}\}\|_{\infty}  = O_{\mathbb{P}} (H  / {N}^{1/2}  ) , \quad
     N  \|{\hat{V}}_{{\hat{\gamma}}}-{\mathbb{E}}\{{\hat{V}}_{{\hat{\gamma}}}\}\|_{\operatorname{op}}   = O_{\mathbb{P}} (H^2  / {N}^{1/2}  ),
$$
which gives sufficient conditions on $H$ to ensure the convergence of ${\hat{V}}_{{\hat{\gamma}}}$ in $L_\infty$ norm and operator norm, respectively.

Finally, equipped with the BEB on the quadratic form $T$ in \eqref{eqn:quad-form} and the conservative variance estimator studied in Theorem \ref{thm:unreplicated-var}, it is immediate to establish Theorem \ref{thm:wald-unreplicated} below for inference, which parallels Theorem \ref{thm:wald-uniform}.

\begin{theorem}[Wald-type inference under unreplicated design]\label{thm:wald-unreplicated}
Consider the unreplicated design with $N_q =1$ for all $q = 1,\ldots, Q$. Assume Conditions \ref{cond:well-conditioned}--\ref{cond:bounded-bgv}.  Then the conclusion of Theorem \ref{thm:wald-uniform} holds if we use the covariance estimator in \eqref{eqn:hV-QU}.
\end{theorem}

\subsection{Non-uniform design}\label{sec::non-uniform-design-inference}

In this section, we consider non-uniform designs in Definition \ref{def:non-uniform-design}. First, we show a BEB on $T$ in \eqref{eqn:quad-form} in Theorem \ref{thm:quad-be-non-uniform} below.

\begin{theorem}[Quadratic form BEB for non-uniform designs]\label{thm:quad-be-non-uniform}
Assume complete randomization that satisfies Definition \ref{def:uniform-design}. Assume Conditions \ref{cond:well-conditioned}, \ref{condition::proper}, \eqref{eqn:small-H} and $N = O(Q)$. There exists a universal constant $C>0$, such that
\begin{align}\label{eqn:quad-be-non-uniform}
    \sup_{t\in{\mathbb{R}}} |{\mathbb{P}}(T\le t)-{\mathbb{P}}(T_0\le t)| \le C \frac{\max_{q\in[Q]}M_N(q)^3 }{\{ \min_{q\in{\mathcal{Q}}_\textsc{S}}S(q,q) \}^{3/2}}\cdot \frac{H^{19/4}}{N^{1/2}}.
\end{align}
\end{theorem}

Theorem \ref{thm:quad-be-non-uniform} assumes $N$ has the same order as $Q$, which is helpful to establish the root $N$ convergence of BEB. We can relax it if we only need the CLT rather than the BEB with .  For the ease of presentation, we omit the general results. A subtle feature of the upper bound in \eqref{eqn:quad-be-non-uniform} is that $\max_{q\in[Q]}M_N(q)$ is the  maximum value of the $M_N(q)$'s over all treatment arms whereas $\min_{q\in{\mathcal{Q}}_\textsc{S}}S(q,q)$ is the minimum value of the $S(q,q)$'s  over treatment arms in ${\mathcal{Q}}_\textsc{S}$ only.

Second, we construct a covariance estimator. It is a combination of the covariance estimators discussed in Sections \ref{sec:var-uniform} and \ref{sec:var-unreplicate}. For the treatment arms with replications, we can calculate sample variances of the potential outcomes based on the observed data. For the treatment arms without replications ${\mathcal{Q}}_\textsc{u}$, we need the grouping strategy in Definition \ref{def::Grouping}. Therefore, we construct a diagonal covariance estimator ${\hat{V}}_{\hat{Y}}$ with the $q$-th diagonal term
\begin{align*}
    {\hat{V}}_{\hat{Y}}(q,q) = \left\{
    \begin{array}{cc}
        \mu_{\langle g \rangle} (Y_q - {\hat{Y}}_{\langle g \rangle})^2, &\quad  q\in{\mathcal{Q}}_{\textsc{u}} \\
        (N_q-1)^{-1}\sum_{Z_i=q}(Y_i-{\hat{Y}}_q)^2, & \quad  q\in{\mathcal{Q}}_{\textsc{r}}\cup{\mathcal{Q}}_{\textsc{l}}.
    \end{array}
    \right.
\end{align*}
In a matrix form, it is equivalent to
\begin{align}\label{eqn:composite-var}
    {\hat{V}}_{\hat{Y}} = \begin{pmatrix}
    {\hat{V}}_{{\hat{Y}},\textsc{u}}& 0& 0,\\
    0& {\hat{V}}_{{\hat{Y}},\textsc{r}}& 0,\\
    0& 0& {\hat{V}}_{{\hat{Y}},\textsc{l}}
    \end{pmatrix},
\end{align}
where ${\hat{V}}_{{\hat{Y}},\textsc{u}}, {\hat{V}}_{{\hat{Y}},\textsc{r}}, {\hat{V}}_{{\hat{Y}},\textsc{l}}$ corresponds to the diagonal covariance estimators for treatment arms ${\mathcal{Q}}_{\textsc{u}}, {\mathcal{Q}}_{\textsc{r}}, {\mathcal{Q}}_{\textsc{l}}$, respectively. Partition the contrast matrix $F$ accordingly as
$$
    F =
    \begin{pmatrix}
    F_{\textsc{s}}\\
    F_{\textsc{l}}
    \end{pmatrix}
    \quad \text{ where }
      F_{\textsc{s}} =
    \begin{pmatrix}
    F_{\textsc{u}}\\
    F_{\textsc{r}}
    \end{pmatrix}.
$$
Construct the final covariance estimator below:
\begin{align}\label{eqn:composite-var-2}
    {\hat{V}}_{\hat{\gamma}} = F^\top {\hat{V}}_{{\hat{Y}}} F = F_{\textsc{u}}^\top {\hat{V}}_{{\hat{Y}},\textsc{u}} F_{\textsc{u}} + F_{\textsc{r}}^\top {\hat{V}}_{{\hat{Y}},\textsc{r}} F_{\textsc{r}} + F_{\textsc{l}}^\top {\hat{V}}_{{\hat{Y}},\textsc{l}} F_{\textsc{l}}.
\end{align}
The decomposition in \eqref{eqn:composite-var-2} allows us to characterize the statistical properties of $ {\hat{V}}_{\hat{Y}} $ by combining the results from Sections \ref{sec:var-uniform} and \ref{sec:var-unreplicate}.

\begin{theorem}[Covariance estimation for non-uniform designs]\label{thm:non-uniform-var}
Consider designs in Definition \ref{def:non-uniform-design} and the covariance estimator in \eqref{eqn:composite-var-2}.
Assume Conditions \ref{cond:moments}, \ref{cond:cond-N}, and \ref{cond:bounded-bgv}.
Assume $\max_{g\in[G]}\mu_{\langle g \rangle}, \Delta$ and $\zeta$ are constants.
\begin{enumerate}[label = (\roman*)]
    \item\label{thm:non-uniform-var-1}

    ${\mathbb{E}}\{{\hat{V}}_{\hat{\gamma}}\}  \succeq V_{\hat{\gamma}}$.

    \item\label{thm:non-uniform-var-2}
    We have

    $$
            \|{\hat{V}}_{{\hat{\gamma}}}-{\mathbb{E}}\{{\hat{V}}_{{\hat{\gamma}}}\}\|^2_{\infty} =
          O_{\mathbb{P}} (   \|F_{\textsc{u}}\|_\infty^4
 |{\mathcal{Q}}_{\textsc{u}}|H^2 +
  \|F_{\textsc{r}}\|_\infty^4 |{\mathcal{Q}}_{\textsc{r}}| H^2  + { \|F_{\textsc{l}}\|_\infty^4 |{\mathcal{Q}}_{\textsc{l}}|^4 N_\textsc{l}^{-3}  H^2} ).
     $$

    \item\label{thm:non-uniform-var-3}
    We have
    $$
            \|{\hat{V}}_{{\hat{\gamma}}}-{\mathbb{E}}\{{\hat{V}}_{{\hat{\gamma}}}\}\|^2_{\operatorname{op}} =
          O_{\mathbb{P}} (   \|F_{\textsc{u}}\|_\infty^4
   |{\mathcal{Q}}_{\textsc{u}}|H^4 +
    \|F_{\textsc{r}}\|_\infty^4 |{\mathcal{Q}}_{\textsc{r}}|  H^4 + { \|F_{\textsc{l}}\|_\infty^4 |{\mathcal{Q}}_{\textsc{l}}|^4 N_\textsc{l}^{-3}  H^4}  )
     $$
\end{enumerate}
\end{theorem}

In Theorem \ref{thm:non-uniform-var}, we assume $\max_{g\in[G]}\mu_{\langle g \rangle}, \Delta$ and $\zeta$ to be constants to simplify the presentation. Without this assumption, we can derive results similar to those in Theorem \ref{thm:unreplicated-var} but relegate finer results to the supplementary material. Theorem \ref{thm:non-uniform-var}\ref{thm:non-uniform-var-1} shows the conservativeness of ${\hat{V}}_{\hat{\gamma}}$ as a direct consequence of Theorems \ref{thm:uniform-var}\ref{thm:uniform-var-1} and  \ref{thm:unreplicated-var}\ref{thm:unreplicated-var-1}. Theorem \ref{thm:non-uniform-var}\ref{thm:non-uniform-var-2} and \ref{thm:non-uniform-var-3} show the stochastic order of the estimation error of ${\hat{V}}_{\hat{\gamma}}$ in $L_\infty$ norm and operator norm, respectively. We only discuss Theorem \ref{thm:non-uniform-var}\ref{thm:non-uniform-var-2} below.
If $\| F \|_{\infty} = O(Q^{-1})$ as in the factorial design in Example \ref{eg::factorial-design}, it reduces to
\begin{align}\label{eqn:V-infty}
   \left\|{\hat{V}}_{{\hat{\gamma}}} - {\mathbb{E}}\{{\hat{V}}_{\hat{\gamma}}\}\right\|_\infty = O_{\mathbb{P}}\left\{ Q^{-2}H(|{\mathcal{Q}}_{\textsc{u}}|^{1/2} + |{\mathcal{Q}}_{\textsc{r}}|^{1/2} + |{\mathcal{Q}}_{\textsc{l}}|^{2}N_{\textsc{l}}^{-3/2})\right\}.
\end{align}
Therefore, if $N$ and $Q$ are of the same order, then $N\|{\hat{V}}_{{\hat{\gamma}}}-{\mathbb{E}}\{{\hat{V}}_{{\hat{\gamma}}}\}\|_{\infty} = O_{\mathbb{P}} ( HN^{-1/2} ) $.
Besides, when one or two of ${\mathcal{Q}}_{\textsc{u}},{\mathcal{Q}}_{\textsc{r}},{\mathcal{Q}}_{\textsc{l}}$ are small or absent, the stochastic orders in Theorem \ref{thm:non-uniform-var} still hold because the large terms in \eqref{eqn:V-infty} will dominate the rest. In particular, if $|{\mathcal{Q}}_\textsc{u}| = |{\mathcal{Q}}_\textsc{r}| = 0$, then
$
    \|{\hat{V}}_{{\hat{\gamma}}} - {\mathbb{E}}\{{\hat{V}}_{\hat{\gamma}}\}\|_\infty = O_{\mathbb{P}} ( HN_{\textsc{l}}^{-3/2}) ,
$
which gives the same rate as Theorem \ref{thm:uniform-var}. If $|{\mathcal{Q}}_\textsc{l}| = 0$, then we should interpret $0\cdot \infty = 0$ in Theorem \ref{thm:non-uniform-var}\ref{thm:non-uniform-var-2} to obtain
$
    \|{\hat{V}}_{{\hat{\gamma}}} - {\mathbb{E}}\{{\hat{V}}_{\hat{\gamma}}\} \|_\infty = O_{\mathbb{P}} ( HQ^{-3/2} )  = O_{\mathbb{P}} ( HN_{\textsc{s}}^{-3/2} ) ,
$
which also agrees with Theorem \ref{thm:uniform-var}.

Finally, the BEB on the quadratic form $T$ and the conservativeness of the covariance estimator ensure Theorem \ref{thm:wald-non-uniform} below for inference.

\begin{theorem}[Wald-type inference under non-uniform designs]\label{thm:wald-non-uniform}
Consider the non-uniform design in Definition \ref{def:non-uniform-design} with \eqref{eqn:small-H} and $N=O(Q)$. Assume Conditions \ref{cond:well-conditioned}--\ref{cond:bounded-bgv}. Then the conclusions of Theorem \ref{thm:wald-uniform}  holds if we use the covariance estimator in \eqref{eqn:composite-var-2}.
\end{theorem}

This concludes our discussion of design-based causal inference with possibly diverging number of treatment levels and varying group sizes across treatment levels.

\section{Simulation}\label{sec:simulation}

In this section, we will evaluate the finite-sample properties of the point estimates and the proposed variance estimator in factorial experiments. We mainly consider non-uniform designs because there have been extensive numerical studies for nearly uniform designs before.

\subsection{Practical implementation}\label{section::simulation-implementation}

For illustration purposes, we focus on conducting inference for the main effects in non-uniform factorial designs.
To do this, we need grouping strategies to implement the proposed variance estimator \eqref{eqn:composite-var}. As we discussed in Section \ref{sec:var-unreplicate}, the structure of factorial designs can provide some practical guidance on the choice of grouping strategy. Besides, our theoretical results in Theorem \eqref{thm:unreplicated-var} also give insights in reducing the conservativeness of the variance estimator. In our simulation, we will compare three variance estimation strategies:
\begin{enumerate}[label = (\roman*)]
    \item \textit{Pairing according to the lexicographical order.} This corresponds to our discussion in Example \ref{exp:pairing}. If arms with similar factor combinations have close means, pairing based on  the lexicographical order can guarantee small between-arm discrepancy in means and reduce the conservativeness.

    Moreover, pairing strategies have another benefit in factorial experiments. We can use a smaller correction factor $\tilde{\mu}_{\langle g\rangle}$ for variance estimation if our goal is to conduct inference marginally (i.e. build confidence intervals on each of $\gamma_h$ separately). The reason is that, while it is hard to control the $\Omega$ matrix in Theorem \ref{thm:unreplicated-var}\ref{thm:unreplicated-var-1} in general, we can control the diagonals of $F_{\textsc{u}}^\top \Omega F_{\textsc{u}}$ because $F_\textsc{u}$ has element $\pm Q^{-1}$. We can get more intuition by noticing that
    $
    \sum_{q'\in{\langle g \rangle_q}, q'\neq q} S(q',q')
    $
    is the core of $\Omega(q,q)$ and that the following algebraic fact holds under pairing:
    \begin{align}\label{eqn:fact-on-Omega}
        \sum_{q\in{\mathcal{Q}}_{\textsc{u}}}\sum_{q'\in{\langle g \rangle_q}, q'\neq q} S(q',q')
        &= \sum_{q\in{\mathcal{Q}}_{\textsc{u}}} S(q,q).
    \end{align}
The identity    \eqref{eqn:fact-on-Omega} enables us to transform the diagonals of $F_{\textsc{u}}^\top \Omega F_{\textsc{u}}$ from source of conservativeness to part of the true variance. Hence it allows us to choose a smaller correction factor:
\begin{align*}
    \tilde{\mu}_{\langle g \rangle} = (1-|{\langle g \rangle}|^{-1})^{-1}(1-3N^{-1})^{-1} = 2(1-3N^{-1})^{-1} ,
\end{align*}
which is approximately one half of the general $\mu_{\langle g \rangle}=4(1-2N^{-1})^{-1}$ given by \eqref{eqn:mu-bg} when $N$ is large.

    \item \textit{Regression-based variance estimation with the target factors as regressors.} Regression-based approach is a commonly used strategy for analyzing factorial experiments \citep{zhao2021regression}.  For non-uniform designs, \cite{zhao2021regression} pointed out ordinary least squares (OLS) with unsaturated model specifications can give biased point estimates and variance estimator. Instead, one should apply weighted least squares (WLS) and the sandwich variance estimation.

    \item \textit{Regression-based variance estimation with the target factors and their high order interactions as regressors.} This strategy differs from strategy (ii) in whether the interactions are included. If all possible two-way interactions of the target factors are specified in the regression model and the true $k$-way ($k\ge 3$) interactions are zero, then this strategy is equivalent to the general factor-based grouping strategy introduced in Example \ref{exp:regression}.
\end{enumerate}

In next section, we will provide more details on implementing the above strategies in simulation.

\subsection{Simulation settings}
We set up a $2^{10}$ non-uniform experiment ($K=10$) following Definition \ref{def:non-uniform-design}, with the basic parameters specified as follows:
\begin{itemize}
    \item Unreplicated arms: $|{\mathcal{Q}}_{\textsc{u}}| = 660$ and $N_q = 1$ for each $q\in{\mathcal{Q}}_{\textsc{u}}$.
    \item Replicated small arms: $|{\mathcal{Q}}_{\textsc{r}}| = 350$ and $N_q = 2$ for each $q\in{\mathcal{Q}}_{\textsc{r}}$.
    \item Large arms: $|{\mathcal{Q}}_{\textsc{l}}| = 14$ and $N_q = 30$ for each $q\in{\mathcal{Q}}_{\textsc{l}}$.
\end{itemize}
The above setup results in a population with $N = 1780$ units. Generate the potential outcomes are independently from a shifted exponential distribution:
\begin{align*}
    Y_i(q) \sim \text{EXP}(\lambda_q) - 1/\lambda_q + \mu_q,
\end{align*}
where $\lambda_q$ are randomly set as $1$ or $2$ with equal probability to induce heteroskedasticity.
We generate two sets of $\mu_q$ and set up two numerical studies, one with small factorial effects and the other with large effects. In both experiments, the main effects for factor $F_k$ with $k=1,4,7,10$ are set as zero. A random subset of two-way of interactions are set as zero as well. All the $k$-way ($k\ge 3$) interactions are zero.
We run simulation for two studies below.
\begin{description}
    \item[Study 1.] \label{study:small} Generated the nonzero main effects and two-way interactions randomly from $ \text{Unif}([-0.5,-0.1]\cup[0.1,0.5])$.
    \item[Study 2.] \label{study:large} Generated the nonzero main effects are randomly from $\text{Unif}([-1,-0.5]\cup[0.5,1])$. Generated the nonzero two-way interactions in the same way as Study 1.
\end{description}

In each study, we focus on estimating the main factorial effects for factor $F_{2l}$ for $l=1,\dots,5$. We apply the point estimates ${\hat{\gamma}}$ in \eqref{eqn:estimates} and compare three variance estimation strategies discussed in Section \ref{section::simulation-implementation} above:
\begin{enumerate}
\item
LEX: We use the grouping strategy based on pairing by the lexicographical order.
\item
$\text{WLS}_0$: We use the sandwich variance estimators based on
WLS with the target factors:
\begin{align*}
    Y \sim F_2 + F_4 + F_6 + F_8 + F_{10}, \text{ with weights } w_i = N_q^{-1}.
\end{align*}
\item
$\text{WLS}_1$: We use the sandwich variance estimators based on   based on WLS with the target factors and their two-way interactions:
\begin{align*}
    Y \sim F_2 + F_4 + F_6 + F_8 + F_{10} + \text{Interaction}_2(F_2, F_4, F_6, F_8, F_{10}), \text{ with weights } w_i = N_q^{-1}.
\end{align*}

\end{enumerate}

\subsection{Simulation results}
We repeat $1000$ times for each study. Figure \ref{fig:simulation} shows the violin plots of the differences between the point estimates and the true parameters. Table \ref{tab:cover-reject} compares the aforementioned variance estimators based on two criteria: coverage rate of $95\%$ confidence intervals and rejection rate of the null that the main effects are zero, which corresponds to the ``Coverage'' column and the ``Rejection'' column, respectively.

Figure \ref{fig:simulation} shows that, even in a highly non-uniform design, the point estimates are centered around the truth and  asymptotic Normality holds as long as the total population $N$ is sufficiently large. Table \ref{tab:cover-reject} shows that the constructed confidence intervals based on all three variance estimators are valid and robust for both small effects and large effects settings. The variance estimator based on pairing is less conservative than the sandwich variance estimator, because the between group variation induced by grouping tends to be smaller with finer groups (see Theorem \ref{thm:unreplicated-var} and the relevant discussion). For the sandwich variance estimator, including more terms into the regression can mitigate the conservativeness. In terms of the rejection rate,  the power of the tests are high for large effects in spite of the conservativeness of the variance estimation. However, if the effects are too small (around two or three times of the true standard deviation of the point estimates), the rejection rate could be negatively impacted.

\begin{figure}[!ht]
    \centering
    \subfigure[Non-uniform experiment with small effects]{
        \includegraphics[width=0.45\textwidth]{BE_PCLT_SMALL_EFFECTS_EXPN.pdf}
        \label{fig:small-effects}
    }
    \subfigure[Non-uniform experiment with large effects]{
        \includegraphics[width=0.45\textwidth]{BE_PCLT_LARGE_EFFECTS_EXPN.pdf}
        \label{fig:large-effects}
    }
    \caption{Violin plots of the differences between the estimators and true parameters for the five
    target effects. Figure \ref{fig:small-effects} corresponds to Study 1 and  Figure \ref{fig:large-effects} corresponds to Study 2, respectively. }
    \label{fig:simulation}
\end{figure}

\begin{table}[!htbp]
\centering
\caption{Coverage and rejection rates based on three variance estimators}
\label{tab:cover-reject}
\begin{tabular}{|c|c|ccc|ccc|}
\hline
\multirow{2}{*}{} & \multirow{2}{*}{Effects} & \multicolumn{3}{c|}{Coverage} & \multicolumn{3}{c|}{Rejection} \\ \cline{3-8}
                         &          & LEX   & WLS-0 & WLS-1 & LEX   & WLS-0 & WLS-1 \\ \hline
\multirow{5}{*}{Study 1} & $F_2$    & 0.963 & 0.977 & 0.973 & 1.000 & 1.000 & 1.000 \\
                         & $F_4$    & 0.968 & 0.980 & 0.976 & 0.026 & 0.015 & 0.018 \\
                         & $F_6$    & 0.974 & 0.985 & 0.982 & 0.665 & 0.600 & 0.622 \\
                         & $F_8$    & 0.974 & 0.984 & 0.982 & 1.000 & 1.000 & 1.000 \\
                         & $F_{10}$ & 0.973 & 0.985 & 0.980 & 0.032 & 0.016 & 0.018 \\ \hline
\multirow{5}{*}{Study 2} & $F_2$    & 0.977 & 0.996 & 0.995 & 1.000 & 1.000 & 1.000 \\
                         & $F_4$    & 0.970 & 0.996 & 0.995 & 0.030 & 0.004 & 0.005 \\
                         & $F_6$    & 0.969 & 0.994 & 0.993 & 1.000 & 1.000 & 1.000 \\
                         & $F_8$    & 0.974 & 0.994 & 0.993 & 1.000 & 1.000 & 1.000 \\
                         & $F_{10}$ & 0.969 & 0.996 & 0.995 & 0.031 & 0.004 & 0.005 \\ \hline
\end{tabular}%
\end{table}

\section{Discussion}\label{sec::discussion}

We focused on scalar outcomes. Results for vector outcomes are also important in both theory and practice. \citet{li2017general} reviewed CLTs and many applications with vector outcomes under the regime of a fixed number of treatment levels and large sample sizes within all treatment levels. We include an extension of the BEB for vector outcomes under a general regime; see Section \ref{sec:vec-outcome} in the supplementary material.

We focused on the theory of linear permutational statistics. Nonlinear permutational statistics also play important roles in studying other problems under the randomization model \citep[e.g.,][]{freedman2008randomization, guo2021generalized}. We will communicate the results in a separate report due to its distinction from the theory in this paper.

Asymptotic results for design-based inference are often criticized  because the population of interest is finite but the asymptotic theory requires a growing sample size. Establishing BEBs is an important theoretical step to characterize the finite-sample performance of the statistics. Alternatively, it is also desirable to derive non-asymptotic concentration inequalities for the estimators under the randomization model. This requires deeper understandings of sampling without replacement and permutational statistics. We leave it to future research.

\bibliographystyle{asa.bst}