The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
37,936 characters
Sharp Minimax Theory for Randomized Experiments
\maketitle
\begin{abstract}
We study minimax-optimal designs and estimators for estimating the sample average treatment effect in finite population randomized experiments, where both design and estimator are unrestricted. For binary potential outcomes, we show this minimax risk is equivalent to the minimax risk $\rho_n^*$ of an estimation problem with $2$ unknown parameters. We leverage this reduction to establish a second-order risk expansion $\rho_n^* = n^{-1} - Cn^{-4/3} + o_n(n^{-4/3})$ for an explicit constant $C$ related to the Airy function. The minimax risk is attained by Bernoulli randomization with a nonlinear shrinkage estimator. Our results show that standard procedures such as complete randomization with difference in means are only minimax optimal up to first order in $n.$ We derive further results on admissibility of these procedures and discuss the practical implications of our results.
\end{abstract}
\noindent\textbf{Keywords:} Minimax risk, experimental design, Bernoulli randomization, nonlinear shrinkage, design-based causal inference, least-favorable prior
\section{Introduction}
\label{sec:introduction}
The fundamental goal of randomized experiments is to estimate the effect of an intervention on a population of interest \cite{fisher1935design}. Two factors are crucial to this enterprise: the choice of experimental design used to randomize the treatment and control, and the choice of estimator used to analyze the data. Naturally, a vast literature has developed on the optimality of experimental designs, and separately, estimators, under various criteria \cite{elfving1952optimum, kiefer1974general, pukelsheim2006optimal, wu1981robustness, li1983minimaxity, hooper1989minimaxity, waite2022minimax, karwa2023admissibility, bai2022optimality, basse2023minimax, kallus2018optimal, kallus2021optimality}. However, a fundamental understanding of optimal \textit{joint} selection over the space of feasible estimators and designs remains incomplete.
An important goal towards closing this gap is to characterize the minimax risk of estimating the sample average treatment effect, optimizing over all choices of estimators and designs. Understanding the minimax risk is important: it serves as a theoretical benchmark of robustness and efficiency, allowing us to distinguish between designs and estimators. This question has been investigated only recently in \cite{kandiros2026design} for finite-population network experiments and in special cases by \cite{bai2023randomize, kallus2018optimal, kallus2021optimality}. In this paper, we provide a precise characterization of this minimax risk under no interference. We determine the sharp constant and second-order term in the minimax risk and show that an optimal choice of design and estimator is given by balanced Bernoulli randomization with a nonlinear shrinkage estimator.
\subsection{Main Results}
Following the potential outcome framework \cite{neyman1923application, rubin1974estimating}, let $Y_i(a), i=1,\dots,n$ denote a fixed binary potential outcome had unit $i$ received treatment $a=0,1$. Let $\tau=\frac{1}{n} \sum_{i=1}^n \left(Y_i(1) - Y_i(0)\right)$ be the sample average treatment effect, and $A_1,\dots,A_n$ be treatment indicators with joint law $\cl{A}$, which we call the \textit{design}. We observe data $Y_i = A_iY_i(1) + (1-A_i)Y_i(0)$ and wish to estimate $\tau$ using an estimator $\hat \tau(A,Y)$. We adopt the terminology in \cite{kandiros2026design} and call a joint choice of design and estimator a \textit{procedure}. For any procedure $(\cl{A},\hat \tau)$, let $R_n(\cl{A},\hat \tau)$ denote its risk $\E[(\hat \tau(A,Y) - \tau)^2]$, where the expectation is taken over the randomness in the design. The choice of risk as mean squared error averaging over the treatment allocation is a standard metric of accuracy, because it succinctly summarizes robustness and efficiency of a given estimator under a fixed design. It is therefore natural to characterize its optimal behavior, in particular the minimax risk
\begin{equation}
\label{eq:minimax_risk_joint_estimator_design}
\inf_{\cl{A}, \hat \tau} \sup_{ \cl{P} \in \set{0,1}^{2n}} R_n(\cl{A},\hat \tau).
\end{equation}
Here, $\cl{P} = \set{Y_i(1),Y_i(0)}_{i=1}^n$ denotes a configuration of potential outcomes and the infimum is taken over all choices of procedure.
We show that \eqref{eq:minimax_risk_joint_estimator_design} admits a characterization as the minimax risk $\rho_n^*$ of a reduced model in which we have nonnegative integer parameters $p,q,r$ which sum to $n$, we observe data $X = p + \Bin(r,1/2)$, and we wish to estimate $(p-q)/n$. For moderate $n$, it is feasible to exactly compute $\rho_n^*$ as well as optimal procedures. A minimax optimal procedure is Bernoulli randomization with equal treatment probability and a nonlinear function applied to the transformed data $\sum_{i=1}^n (A_iY_i + (1-A_i)(1-Y_i)).$ The nonlinearity corresponds to shrinkage towards zero under a least-favorable prior. These results are discussed in Section \ref{sec:reduction_to_simplified_model}, which also establishes similar results for bounded potential outcomes whether discrete or continuous.
On the theoretical front, we establish in Section \ref{sec:second_order_risk} a second-order expansion of the minimax risk
\begin{equation}
\label{eq:second_order_risk_minimax_intro}
\rho_n^* = \frac{1}{n} - \frac{C_A}{n^{4/3}} + o_n(n^{-4/3}).
\end{equation}
The constant $C_A$ is related to the Airy function, one of the linearly independent solutions of the differential equation $d^2y/dx^2 - yx = 0.$ The asymptotic expansion to second-order is evocative of similar results in the minimax estimation of bounded normal means \cite{bickel1981minimax} and other settings \cite{johnstone1992minimax}. The proof strategy for \eqref{eq:second_order_risk_minimax_intro} is similar to the line of work by Levit \cite{levit1981asymptotic,levit1983minimax,levit1986second}, which exposes a connection between the second-order term in a minimax risk expansion and the solution of a problem-specific differential equation, in this case the Airy equation. To our knowledge, the appearance of the $n^{-4/3}$ term and connection to the Airy function is nonstandard and new in this literature.
Our results provide novel perspectives on standard procedures, carefully discussed in Section \ref{sec:comparison_to_standard_methods}. Common procedures in the literature typically use complete randomization (\textsf{CRE}) or Bernoulli randomization (\textsf{BRE}). These designs appear in the basic textbooks in the field \cite{imbens2015causal, ding2024first} and have a long, rich tradition \cite{fisher1992arrangement,neyman1923application,kempthorne1955randomization}. For example, Imbens and Rubin name these among the four classical experimental designs \cite{imbens2015causal}. These designs are often paired with the difference-in-means (\textsf{DIM}) or the Horvitz-Thompson (\textsf{HT}) estimators. Much work in the design-based causal inference literature focuses on these estimators \cite{aronow2013class,aronow2014sharp,lin2013agnostic, athey2017econometrics, imbens2015causal, chattopadhyay2024neymanian, li2017general, harshaw2021algorithmic}. In light of \eqref{eq:second_order_risk_minimax_intro}, we find that the procedures (\textsf{CRE}, \textsf{DIM}), (\textsf{BRE}, \textsf{DIM}) are only minimax-optimal to first order. Figure \ref{fig:minimax_risk_curves} displays the maximum risks of these procedures alongside $\rho_n^*$. For small to moderately sized $n$, the maximum risk of these procedures is substantively larger than $\rho_n^*$, which also serves as a lower bound for the maximum risk. For example, at $n = 100$, the maximum risk for (\textsf{CRE}, \textsf{DIM}) is about 40\% larger than the optimum. Randomized controlled trials in practice are often of this size \cite{marshall2021state}. Complementing our analysis on minimaxity, we show that (\textsf{BRE}, \textsf{Opt}) is admissible. On the other hand, (\textsf{CRE}, \textsf{DIM}) is also admissible and dominates two other commonly used procedures, (\textsf{BRE}, \textsf{DIM}), (\textsf{BRE}, \textsf{HT}). Neither (\textsf{BRE}, \textsf{Opt}) nor (\textsf{CRE}, \textsf{DIM}) dominates the other.
\begin{figure}[H]
\centering
\includegraphics[width=0.8\linewidth]{figures/scaled_minimax_risk_comparison.png}
\caption{Maximum Risks of Various Designs and Estimators, multiplied by $n$. The maximum risk of \textsf{BRE} \ with the Horvitz-Thompson estimator is $4/n$ and not plotted here.}
\label{fig:minimax_risk_curves}
\end{figure}
\subsection{Related Work}
Most closely related to this work are previous studies by Kallus \cite{kallus2018optimal,kallus2021optimality}, Bai \cite{bai2023randomize}, and Kandiros et al. \cite{kandiros2026design}. These works have a slightly richer model of nature than the present paper, allowing for covariates, distributions on potential outcomes, or violations of SUTVA. Thus, these works are complementary to our paper. In \cite{kallus2018optimal,kallus2021optimality}, Kallus solves for the minimax optimal design, fixing the estimator to be difference in means. Complete randomization is shown to be minimax optimal under bounded norm assumptions on the potential outcome matrix. Under restrictions on the relationship between the potential outcomes and covariates, other designs are found to be minimax optimal. \cite{bai2023randomize} generalizes some of the previous results. Under a superpopulation framework, Bai proves a general design minimaxity theorem for estimating arbitrary contrasts of potential outcome means, assuming a class of potential outcome distributions invariant under some group of permutations. As a corollary of this general theorem, $(\textsf{CRE},\textsf{DIM})$ is found to be the minimax optimal procedure restricted to estimators which are \textit{linear} in the observed data $Y$.
Kandiros et al. \cite{kandiros2026design} tackle our same question of minimax optimality over joint choices of design and estimator for network interference experiments. They introduce the terminology \textit{procedure} to refer to such a joint choice, which we will follow throughout this paper. \cite{kandiros2026design} characterizes optimal minimax rates for estimating the global average treatment effect and other estimands, with upper and lower bounds that depend on properties of the interference graph. Specializing to the SUTVA case gives a minimax rate of $\Theta(1/n)$. The present paper sharpens the SUTVA analysis, precisely characterizing the constant as well as the second-order term. In the interference literature, \cite{karwa2023admissibility} provides another example of decision-theoretic optimality results on standard causal estimators.
The problem of determining a minimax procedure has also been considered in survey sampling \cite{stenger1979minimax, stenger1988asymptotic, rinott2009some}. Proposition 17 of the survey by Rinott \cite{rinott2009some} provides the answer: Among procedures sampling $n$ units from a population of size $N$, the minimax optimal procedure is simple random sampling with an affine shrinkage rule introduced in \cite{hodges1982minimax}. The minimax procedure and risk admit explicit formulas, from which it can be shown that
\begin{equation}
\frac{a_\lambda}{N} - \frac{b_\lambda}{N^{3/2}} + o(N^{-3/2}),
\end{equation}
for a finite population with elements in $[0,1]$ and $n = \floor{\lambda N}, \lambda \in (0,1)$, with some constants $a_\lambda,b_{\lambda}.$ This result serves as an interesting point of comparison to the randomized experiment setting. The second-order term decays faster than in the randomized experiment setting, and
affine shrinkage is optimal instead of nonlinear shrinkage.
The recent paper
\cite{aronow2026minimax} extends these results, deriving the minimax optimal choice of procedure under unbiasedness restrictions and non-symmetric parameter spaces. Therein, they show the centered version of the Horvitz-Thompson estimator \cite{aronow2013class} is minimax optimal. This estimator also plays a role in the randomized experiment setting. We also note that results separately optimizing the sampling design or estimator abound in this literature \cite{godambe1965admissibility, hodges1982minimax, cheng1983minimax, cheng1987optimality, bickel2011minimax, joshi1979best, scott1975minimax}.
\section{Reduction to a Simplified Model}
\label{sec:reduction_to_simplified_model}
In this section, we describe the reduced two-parameter model with minimax risk $\rho_n^*$ and analyze its properties. The connection to the original procedure problem is made in the next subsection. Consider first a statistical model with unknown parameters $\Delta_1,\dots,\Delta_n \in \set{-1,0,1}$ and data $(S_1,\dots,S_n)$ given by
\begin{equation}
\label{eq:S_definition}
\begin{cases}
S_i = 1 & \text{ if } \Delta_i = 1 \\
S_i = 0 & \text{ if } \Delta_i = -1 \\
S_i \sim \textsf{Ber}(1/2) & \text{ if } \Delta_i = 0,
\end{cases}
\end{equation}
all mutually independent. $\Delta_i$ will be interpreted as the individual treatment effects $Y_i(1) - Y_i(0)$ in the context of the original problem. Let $p$ denote the number of indices $i$ with $\Delta_i = 1,$ $q$ the number of indices with $\Delta_i = -1$ and $r$ the remainder. We are interested in estimating
\[
\tau(\Delta) := \frac{1}{n}\sum_{i=1}^n \Delta_i = \frac{p - q}{n}.
\]
Let $r_n^*$ denote the minimax risk for estimating $\tau(\Delta)$ given data $(S_1,\dots,S_n)$. By an invariance argument, it suffices to keep track of the parameters $p,q,r$ and to use estimators which are functions of $\sum_{i=1}^n S_i$. Therefore, it suffices to study the minimax risk of the even simpler model with parameters $(p,q,r)$ in the index set
\begin{equation}
\cl{I}_n := \set{(p,q,r) \in \bb{Z}_{\geq 0}^3: p+q+r = n},
\end{equation}
where we observe data $X \sim p + \Bin(r,1/2)$, with estimand $(p-q)/n$. Let $\rho_n^*$ be the minimax risk of this binomial experiment, which we will refer to as the \textit{reduced model} throughout. Then:
\begin{theorem}
\label{thm:simplified_model_minimax_risk}
\begin{equation}
\label{eq:simplified_minimax_risk}
r_n^* = \rho_n^* = \inf_f \sup_{\cl{I}_n} \E\left(f(p + \textsf{Bin}(r,1/2)) - \frac{p-q}{n} \right)^2.
\end{equation}
Moreover, there is a unique $f^*_n$ minimizing \eqref{eq:simplified_minimax_risk} and it is the posterior mean of $(p-q)/n$ under a least-favorable prior $\pi_n^*$ on $\cl{I}_n$. As a consequence, $f^*_n$ is admissible and $\rho_n^*$ is the Bayes risk of $\pi_n^*.$
\end{theorem}
Although the minimax estimator $f_n^*$ is unique, the least favorable prior need not be unique. Computing the minimax optimal estimator and least-favorable prior in the reduced model can be done using convex optimization for moderate values of $n$ below a few hundred. In contrast, directly attempting to solve \eqref{eq:minimax_risk_joint_estimator_design} is computationally prohibitive. Figure \ref{fig:numerical_solutions_minimax} shows the numerical solutions $f_n^*$ to \eqref{eq:simplified_minimax_risk}, rewritten as functions of the natural unbiased estimator $(2\sum_{i=1}^n S_i - n)/n$. It is clear from the right-hand panel of the figure that the optimal minimax estimators nonlinearly shrink the unbiased estimator $(2\sum_{i=1}^n S_i - n)/n$ towards zero.
\begin{figure}
\centering
\includegraphics[width=\linewidth]{figures/numerical_minimax_estimators.pdf}
\caption{Numerical solutions of the minimax estimator for various values of $n$, displayed as functions of $(2\sum_{i=1}^n S_i - n)/n$. In particular, $g_n^*(x) := f_n^*\left( n(x + 1)/2\right)$ is plotted. The right panel shows deviation from the identity and emphasizes the nonlinear shrinkage.}
\label{fig:numerical_solutions_minimax}
\end{figure}
The equivalence between the reduced model and randomized experiments is formally made by treating the transformed data $A_iY_i + (1 - A_i)(1 - Y_i)$ as the observations $S_i$ and interpreting $\Delta_i$ as the ITE $Y_i(1) - Y_i(0)$. These quantities $S$ have been considered before in \cite{jaskowski2012uplift}, where they are called \textit{class-variable transformations} and utilized for heterogeneous treatment effect targeting.
\begin{theorem}
\label{thm:minimax_risk_reduction}
Let $f_n^*$ be the function minimizing \eqref{eq:simplified_minimax_risk} in the reduced model. Then
\begin{equation}
\label{eq:minimax_risk_is_equal_to_simplified_minimax_risk}
\inf_{\cl{A}, \hat \tau} \sup_{\cl{P} \in \set{0,1}^{2n}} R_n(\cl{A},\hat \tau) = \rho_n^*.
\end{equation}
A minimax procedure achieving the upper bound is given by Bernoulli$(1/2)$ randomization with corresponding estimator $\hat{\tau}_{\textsf{Opt}}=f^*_n(\sum_{i=1}^n S_i)$. Furthermore, this procedure is admissible.
\end{theorem}
We will refer to this minimax admissible procedure as (\textsf{BRE}, \textsf{Opt}). A proof sketch of Eq. \eqref{eq:minimax_risk_is_equal_to_simplified_minimax_risk} is instructive. The idea for the lower bound is to consider a particular prior on the potential outcomes $Y_i(1),Y_i(0)$, which exactly mirrors the reduced model. For any fixed vector $\Delta \in \set{-1,0,1}^n$, generate $(Y_i(1),Y_i(0))$ according to
\begin{equation}
\label{eq:bayes_model}
(Y_i(1),Y_i(0)) =
\begin{cases}
(1,0) & \text{ if } \Delta_i = 1 \\
(0,1) & \text{ if } \Delta_i = -1 \\
(0,0) \text{ or } (1,1) \text{ w.p. } 1/2 & \text{ if } \Delta_i = 0.
\end{cases}
\end{equation}
The crucial property is that the transformed quantities $S_i = A_iY_i(1) + (1 - A_i)(1-Y_i(0))$ are independent of the design $A$, under this model. There is a bijection between $(A,Y)$ and $(A,S)$, and so estimators which depend only on $S$ may be considered. The optimal estimator is then the posterior mean of $\tau,$ which equals $\frac{1}{n}\sum_{i=1}^n \Delta_i$, given the data $S$. This is the Bayes risk of the reduced model under a particular prior on $\Delta$. Choosing a least-favorable prior and applying a finite minimax theorem shows the lower bound. The upper bound is achieved simply by taking the $\textsf{BRE}(1/2)$ design, under which the quantities $S_i$ have exactly the law in \eqref{eq:S_definition}.
Although the setting of binary outcomes is already highly relevant in causal inference \cite{ding2024first, imbens2015causal}, the results can be extended further to bounded potential outcomes. The next result characterizes the minimax risk when the potential outcomes are known to lie in an interval $[L,U].$
\begin{theorem}
\label{thm:minimax_risk_reduction_bounded_case}
For any finite $U > L$, we have
\begin{equation}
\inf_{\cl{A}, \hat \tau} \sup_{\cl{P} \in [L,U]^{2n}} R_n(\cl{A},\hat \tau) = (U - L)^2r_n^*.
\end{equation}
A minimax optimal procedure is $\textsf{BRE}(1/2)$ with the estimator
$\E\left[f^*_n\left(\sum_{i=1}^n B_i\right) \mid S \right], B_i \stackrel{ind}{\sim} \Ber(S_i).$
\end{theorem}
If the outcomes are unbounded, the theorem immediately implies that the worst-case risk is infinite.
\section{Analysis of Minimax Risk}
\label{sec:second_order_risk}
Figure \ref{fig:minimax_risk_curves} shows the curious fact that $n\rho_n^*$ increases rather slowly to its limiting value of $1$. The next theorem, our main result, provides an explanation by expanding the minimax risk up to second order. It also provides insight into the form of the optimal minimax estimator and the associated least-favorable prior.
\begin{theorem}
\label{thm:second_order_minimax_risk_behavior}
Let $a_1'$ be the largest negative zero of the derivative of the Airy function. The minimax risk satisfies
\[
\rho_n^* = \frac{1}{n} - \frac{C_A}{n^{4/3}} + o_n(n^{-4/3}),
\]
with $C_A = -4^{1/3}a_1' = 1.617\dots$.
\end{theorem}
To give some intuition for the proof, we first reparametrize the model of Section \ref{sec:reduction_to_simplified_model} by considering $2\sum_i S_i - n = \sum_{i=1}^r \e_i + (p-q)$ with $\e_i$ i.i.d. Rademacher $1/2$ random variables. Let $\theta := (p-q)$. Then $\rho_n^*$ is equivalent to the minimax risk of estimating $\theta/n$ on the basis of the observation $X = \left(\theta + \sum_{i=1}^r \e_i\right)$ with unknown parameters $(\theta,r)$ in the space
\begin{equation}
\Theta_n := \set{(\theta,r): 0 \leq r \leq n, |\theta| + r \leq n, \theta + r \equiv n \mod 2}.
\end{equation}
For each $(\theta,r) \in \Theta_n$ there exists $p,q,r$ such that $p - q = \theta$ and $p,q,r \geq 0 , p + q + r = n,$ so this is indeed an equivalent parametrization.
In this model, it is not difficult to show that $\rho_n^* = (1 + o(1))/n$ with the natural unbiased estimator $X/n$ being first-order minimax with maximum risk $1/n$. Therefore, we expect estimators achieving minimax risk up to second order to be small perturbations of the estimator $X/n$. Consider the candidate estimator
\begin{equation}
\label{eq:second_order_minimax_estimator_h}
\delta_{n,h}(X) := \frac{X}{n} - n^{-\alpha} h(X/n^\alpha),
\end{equation}
for some function $h$ and $\alpha \geq 0$. To heuristically choose $h$ and $\alpha$, suppose $h$ is sufficiently well-behaved so that $h'$ and $h^2$ are Lipschitz. Then the risk $\delta_{n,h}(X)$ can be expanded by using Lemma \ref{lemma:stein_rademacher_sums}, a version of Stein's identity for i.i.d. sums of Rademacher random variables. Let $U_r = \sum_{i=1}^r \e_i$, $\tilde{\theta} = \theta/n^\alpha$. We have
\begin{align*}
\E[(\delta_{n,h}(X) - \theta/n)^2] & = \E\left[\left(\frac{U_r}{n} - n^{-\alpha} h(X/n^{\alpha})\right)^2 \right] \\
& = \frac{r}{n^2} - 2n^{-\alpha - 1}\E[U_r h(\tilde \theta + U_r/n^\alpha)] + n^{-2\alpha}\E[h(\tilde \theta + U_r/n^\alpha)^2] \\
& \approx \frac{r}{n^2} - 2rn^{-2\alpha - 1} \E[h'(\tilde{\theta} + U_{r}/n^{\alpha})] + n^{-2\alpha}h(\tilde \theta)^2 + O(n^{-3\alpha}\sqrt{r}) \\
& = \frac{r}{n^2} - 2rn^{-2\alpha - 1} h'(\tilde \theta)+ n^{-2\alpha}h(\tilde \theta)^2 + O(n^{-3\alpha}\sqrt{r}).
\end{align*}
It is intuitively clear that the hardest parameters $(\theta,r)$ are those on the boundary $r = n - |\theta|.$ Thus, ignoring lower order terms, the worst-case risk is given by
\begin{equation}
\sup_{\tilde \theta} \frac{1}{n} - \frac{|\tilde \theta|}{n^{2-\alpha}} - \frac{2h'(\tilde \theta)}{n^{2\alpha}} + \frac{h(\tilde \theta)^2}{n^{2\alpha}}.
\end{equation}
The sharpest upper bound on the minimax risk is then obtained by solving the following problem:
\[
\frac{1}{n} + \inf_h \sup_{x} \left\lbrace \frac{h(x)^2 - 2h'(x)}{n^{2\alpha}} - \frac{|x|}{n^{2-\alpha}} \right \rbrace.
\]
For the latter term to remain second-order, we require $\alpha \in (1/2,1).$ For $\alpha < 2/3$, the second-order term in the minimax problem is $\inf_h \sup_x (h(x)^2 - 2h'(x))/n^{2\alpha}$. It can be shown that the solution is trivial and given by $h = 0.$ If $\alpha > 2/3$, the second-order term is $-|x|/n^{2-\alpha}$ and the perturbation is irrelevant. Therefore, a second-order risk improvement is obtained only when $\alpha = 2/3$, which yields
\begin{equation}
\frac{1}{n} + \inf_h \sup_{x} \left\lbrace \frac{h(x)^2 - 2h'(x) - |x|}{n^{4/3}} \right \rbrace,
\end{equation}
and the variational problem of choosing $h$ to minimize $M(h) := \sup_{x \in \bR} h(x)^2 - 2h'(x) - |x|.$ By integrating against smooth test functions $\phi(x)^2$ and completing the square, it can be shown that
\begin{equation}
\label{eq:variational_problem}
M(h) \geq - \inf_{\norm{\phi}_2 = 1} \int_{\bR} 4\phi'(x)^2 + |x|\phi(x)^2 dx
\end{equation}
with equality if and only if $h(x) = -2\phi_A'(x)/\phi_A(x)$, where $\phi_A$ is the optimizer of the right-hand side. The quantity $\cl{E}(\phi) = \int_\bR 4\phi'(x)^2 + |x|\phi(x)^2$ is the quadratic form of the second-order differential operator $-4\frac{d^2}{dx^2} + |x|$. By the Rayleigh quotient characterization, the value of the variational problem $\inf_\phi \cl{E}(\phi)$ is the smallest eigenvalue of the differential operator. This turns out to be the constant $C_A$ in the statement of Theorem \ref{thm:second_order_minimax_risk_behavior}, so $M(h) \geq - C_A.$ The corresponding eigenfunction $\phi_A$ satisfies $-4\phi_A''(x) + |x|\phi_A(x) = C_A \phi_A(x),$ from which it can be shown that
\begin{equation}
\phi_A(t) \propto \Ai\left( \frac{|t| - C_A}{4^{1/3}}\right).
\end{equation}
Thus, taking $h = -2\phi_A(x)'/\phi_A(x)$ in \eqref{eq:second_order_minimax_estimator_h} yields the candidate second-order minimax estimator. The lower bound establishing Theorem \ref{thm:second_order_minimax_risk_behavior} is obtained by considering priors that are discrete approximations of $\phi_A(x)^2$.
This tight interplay between second-order differential equations and variational problems appears throughout the literature on second-order minimax risk expansions. Essentially the same strategy was exploited for bounded normal mean estimation, pioneered primarily by a sequence of papers by Levit \cite{levit1981asymptotic, levit1983minimax, levit1986second, koranyi2002asymptotically} and also Bickel \cite{bickel1981minimax}. In this literature, the Gaussian structure is exploited by using Brown's identity to relate Bayes risk to prior information. In our proof, we instead use Van Trees' inequality, a softer tool. A similar interplay happens in \cite{johnstone1992minimax, Aslan2006} outside of the Gaussian setting.
The appearance of the Airy function and the slower $n^{-4/3}$ second-order term seems to be novel. The primary reason underlying this phenomenon is that the maximum variance $\Var(X/n) = (n - |\theta|)/n^2 = \frac{1}{n} - \frac{|\theta|}{n^2}$ at a fixed $\theta$ decays linearly in $|\theta|$. In estimation problems with bounded parameter spaces like the bounded normal means problem, the variance is typically constant.
\section{Comparison to Standard Procedures}
\label{sec:comparison_to_standard_methods}
In practice, the most widely used designs are the completely randomized experiment and Bernoulli randomized experiment \cite{ding2024first, imbens2015causal}. Recall that the \textsf{CRE} \ selects $n_1$ of the $n$ units uniformly at random to be treated, and Bernoulli randomization (\textsf{BRE}) assigns treatment via indicators $A_i$ drawn i.i.d. $\Ber(p)$. We will focus on the \textsf{CRE} \ where $n/2$ units are treated, and on \textsf{BRE} s with treatment probability $1/2$. For simplicity throughout this section, we will assume $n$ is even. Two commonly used estimators are the difference-in-means and Horvitz-Thompson estimators
\begin{align}
\label{eq:dim}
\hat \tau_{\textsf{DIM}}(A,Y) & := \frac{\sum_{i=1}^n A_i Y_i}{\sum_{i=1}^n A_i} - \frac{\sum_{i=1}^n (1-A_i) Y_i}{\sum_{i=1}^n (1-A_i)}, \\
\label{eq:horvitz_thompson}
\hat \tau_{\text{\textsf{HT}}}(A,Y) & := \frac{\sum_{i=1}^n A_i Y_i}{n/2} - \frac{\sum_{i=1}^n (1-A_i) Y_i}{n/2}.
\end{align}
In addition, we will consider the centered Horvitz-Thompson estimator \textsf{cHT}, which has been studied in \cite{aronow2013class, aronow2026minimax}:
\begin{equation}
\label{eq:centered_horvitz_thompson}
\hat \tau_{\text{\textsf{cHT}}}(A,Y) := \frac{\sum_{i=1}^n A_i (Y_i - 1/2)}{n/2} - \frac{\sum_{i=1}^n (1-A_i) (Y_i-1/2)}{n/2}.
\end{equation}
$\hat \tau_{\text{\textsf{cHT}}}$ is a special case of the family of estimators introduced in \cite{aronow2013class}. It appears in essentially the same form in the survey sampling literature \cite{aronow2026minimax}, where their form of $\hat \tau_{\text{\textsf{cHT}}}$ is shown minimax optimal over design-unbiased estimators under independent sampling. The estimator arises naturally in our problem: a quick derivation shows that $\hat \tau_{\textsf{cHT}}(A,Y) = \frac{1}{n}\left(2\sum_{i=1}^n S_i - n\right)$, which is exactly the natural unbiased estimator in the reduced model of Section \ref{sec:reduction_to_simplified_model}.
For the \textsf{CRE}, $\textsf{cHT}$ is equal to both $\textsf{DIM}$ and $\textsf{HT}$; in \textsf{BRE}, they are distinct. For (\textsf{BRE}, \textsf{DIM}), we take the convention $\hat \tau_{\textsf{DIM}} = 0$ if $\sum_i A_i \in \set{0,n}.$ Note that in the Bernoulli randomized experiment, \textsf{HT} \ is known to be asymptotically less efficient than \textsf{DIM}, but we retain it nonetheless. We will therefore consider four procedures: (\textsf{CRE}, \textsf{DIM}), (\textsf{BRE}, \textsf{DIM}), (\textsf{BRE}, \textsf{cHT}) and (\textsf{BRE}, \textsf{HT}). By computing their maximum risks, we can show that the three procedures $(\textsf{CRE},\textsf{DIM}), (\textsf{BRE},\textsf{DIM}), (\textsf{BRE},\textsf{cHT})$ are only minimax to first order, while $(\textsf{BRE},\textsf{HT})$ is not.
\begin{proposition}[Maximum Risk]
\label{prop:maximum_risks}
Let $n$ be even. Then the following are true.
\begin{align*}
\sup_{\cl{P} \in \set{0,1}^{2n}} R_n(\textsf{CRE},\textsf{DIM}) & = \frac{1}{n-1} \qquad
&\sup_{\cl{P} \in \set{0,1}^{2n}} R_n(\textsf{BRE},\textsf{HT}) = \frac{4}{n} \\
\sup_{\cl{P} \in \set{0,1}^{2n}} R_n(\textsf{BRE},\textsf{DIM}) & = \frac{1}{n} + \frac{2}{n^2} + O(n^{-3}) \qquad
&\sup_{\cl{P} \in \set{0,1}^{2n}} R_n(\textsf{BRE},\textsf{cHT}) = \frac{1}{n}.
\end{align*}
\end{proposition}
The proof of Proposition \ref{prop:maximum_risks} also exposes the worst-case potential outcome configurations for each procedure. For the procedures using $\textsf{DIM}$ with $n$ even, the worst-case is $n/2$ units with $(Y_i(1),Y_i(0)) = (1,1)$ and $n/2$ units with $(0,0)$. These worst-case configurations have zero SATE and are unique up to relabeling of units. For $(\textsf{BRE},\textsf{HT})$, the worst-case is $Y_i(1) = Y_i(0) = 1$ for all $i;$ for $(\textsf{BRE},\textsf{cHT}),$ every configuration for which $Y_i(1) - Y_i(0) = 0, \forall i$ achieves the maximum risk. For $(\textsf{BRE},\textsf{Opt})$, there are multiple worst-case configurations. Recall that the minimax optimal procedure is a function of $\sum_{i=1}^n S_i$, and under \textsf{BRE}, $\sum_{i=1}^n S_i \sim p + \Bin(r,1/2)$, with $S_i = Y_iA_i + (1-Y_i)(1-A_i)$ and $p,q,r$ being the number of units with ITE equal to $1,-1,0$ respectively. Then any potential outcome configuration with the same $(p,q,r)$ values shares the same risk. Figure \ref{fig:opt_worst_case_configurations} displays some numerically-computed worst-case values of $(p,q,r)$ for various values of $n$. Note there are worst-case configurations with SATE equal to $-1,0,1$ and various values in between. All configurations appear to satisfy $p = 0$ or $q = 0$ and the set is symmetric upon swapping $p,q$.\\
\begin{figure}
\centering
\includegraphics[width=\linewidth]{figures/worst_case_configurations_n10_n20_n50.pdf}
\caption{Blue dots represent numerically computed worst-case configurations for $(\textsf{BRE},\textsf{Opt})$, up to a tolerance of $10^{-7}.$ The plots represent the two-dimensional simplex $p+q+r = n$. Some explicit labels are provided on the left plot.}
\label{fig:opt_worst_case_configurations}
\end{figure}
Despite its aesthetic and theoretical appeal, minimaxity may be too conservative a criterion to have practical implications -- a practitioner may suspect these worst-case configurations are too unrealistic. A useful complement is to consider admissibility and dominance as another lens. Recall that a procedure $(\cl{A},\hat\tau)$ \textit{dominates} another procedure $(\cl{B},\hat\theta)$ if
\[
\E_{\cl{A}}\left[\left(\hat\tau(A,Y) - \tau \right)^2\right] \leq \E_{\cl{B}}\left[\left(\hat\theta(A,Y) - \tau \right)^2\right]
\]
with strict inequality at one configuration $\cl{P} \in \set{0,1}^{2n}.$ A procedure is \textit{admissible} if no other procedure dominates it.
\begin{theorem}[Admissibility of procedures]
\label{thm:admissibility_of_procedures}
For all even $n \geq 4$, (\textsf{CRE}, \textsf{DIM}) and (\textsf{BRE}, \textsf{cHT}) are admissible. Furthermore (\textsf{CRE}, \textsf{DIM}) dominates (\textsf{BRE}, \textsf{DIM}), (\textsf{BRE}, \textsf{HT}). Therefore, they are inadmissible.
\end{theorem}
The domination of (\textsf{BRE}, \textsf{DIM}), (\textsf{BRE}, \textsf{HT}), all else equal, suggests we should prefer (\textsf{BRE}, \textsf{Opt}), (\textsf{BRE}, \textsf{cHT}), or (\textsf{CRE}, \textsf{DIM}). Although (\textsf{BRE}, \textsf{Opt}) demonstrates a substantial improvement in maximum risk per Figure \ref{fig:minimax_risk_curves}, a clear practical recommendation of which admissible procedure to use requires further research. We provide some numerical comparisons of the risk of these 3 procedures across various subsets of the potential outcome space $\set{0,1}^{2n}$.
The first simulation determines the fraction of the configuration space $\set{0,1}^{2n}$ for which each procedure has the smallest risk amongst all competitors. Two measures are used -- the standard uniform measure on the hypercube, and the measure under which count vectors $(n_{11},n_{10},n_{01},n_{00})$ are weighted equally, where $n_{ab} = \#\set{i: Y_i(1)=a,Y_i(0) = b}$. Figure \ref{fig:admissible_procedure_comparison_full_space} shows the results. It preempts the criticism that (\textsf{BRE}, \textsf{Opt}) might only reduce risk at its worst-case configurations, while being uncompetitive everywhere else. Under the uniform prior as $n$ increases, (\textsf{BRE}, \textsf{Opt}) has smaller risk than the other procedures on nearly every configuration. Under the count-vector weighting scheme (\textsf{CRE}, \textsf{DIM}) appears to be the best performer, although (\textsf{BRE}, \textsf{Opt}) is still useful. One reason explaining the phenomenon in the left panel of Figure \ref{fig:admissible_procedure_comparison_full_space} is simply that (\textsf{BRE}, \textsf{Opt}) uses a shrinkage estimator and performs well when the SATE is small. Figure \ref{fig:fixed_tau_admissible_procedure_comparison} corroborates this, showing the same comparisons on the subsets of the configuration space with fixed $\tau.$ As $|\tau|$ increases, (\textsf{BRE}, \textsf{Opt}) is less competitive than the other procedures. Note that these plots do not show the difference in risk between the three methods.
\begin{figure}
\centering
\includegraphics[width=\linewidth]{figures/admissible_procedure_comparison.pdf}
\caption{Each line shows the measure of configuration space $\set{0,1}^{2n}$ on which a fixed procedure has strictly smaller risk than its competitors, among the three admissible procedures for even $n$ studied in this paper. Ties are ignored. Two measures are used as in the text.}
\label{fig:admissible_procedure_comparison_full_space}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=\linewidth]{figures/conditional_parameter_space_win_volume.pdf}
\caption{Same as Figure \ref{fig:admissible_procedure_comparison_full_space}, with configuration space restricted to those with fixed $\tau$. Only values of $n$ for which $n\tau \in \bb{Z}$ and which are multiples of $10$ are plotted.}
\label{fig:fixed_tau_admissible_procedure_comparison}
\end{figure}
\section{Discussion}
\label{sec:discussion}
In this paper, we offer a sharp description of the minimax risk of estimating the sample average treatment effect with bounded potential outcomes, minimizing over joint choices of design and estimator. We show the minimax risk is exactly that of a simplified experiment with parameters $p,q,r \geq 0$ such that $p+q+r = n$, we observe data $X \sim p + \Bin(r,1/2)$, and wish to estimate $(p - q)/n$. We derive a second-order expansion for the minimax risk of this simplified experiment, revealing the slow second-order rate $n^{-4/3}$ and an interesting connection to the Airy function.
Although our results are intended to be a theoretical benchmark, the procedure $(\textsf{BRE},\textsf{Opt})$ may potentially be developed into a practical method for the analysis of small randomized controlled trials. The sample size $n$ should be in some intermediate range (perhaps $n \leq 100$) where computing the exact rule is not computationally prohibitive and the nonlinear rule offers a substantial reduction in maximum risk over $(\textsf{CRE},\textsf{DIM}), (\textsf{BRE},\textsf{cHT})$. Broad surveys of medical RCTs suggest that the number of participants recruited frequently falls within this range \cite{totton2023review, chang2022survey, marshall2021state}. Our numerical results in Figures \ref{fig:admissible_procedure_comparison_full_space}, \ref{fig:fixed_tau_admissible_procedure_comparison} also suggest $(\textsf{BRE},\textsf{Opt})$ is a useful procedure for small effect sizes, which are important in practice. Also, the ability to perform inference with the minimax optimal procedure would be an important problem to address, as well as the potential to incorporate measured baseline covariates to further enhance efficiency and power. For inference, a randomization test might be well-suited for $(\textsf{BRE},\textsf{Opt})$. The test may be inverted to give confidence intervals by considering weak null hypotheses as in \cite{wu2021randomization}.
We put forward (\textsf{BRE}, \textsf{Opt}) as a provably minimax optimal procedure deserving more research to determine its practical applicability. In the meantime, our results justify the continued use of (\textsf{CRE}, \textsf{DIM}) and (\textsf{BRE}, \textsf{cHT}), both of which are first-order minimax optimal, admissible over all procedures, unbiased, and have a well-developed inferential theory.
\section*{Acknowledgements}
TS thanks Lihua Lei, Harrison Li, and Maggie Wang for helpful comments. This work was supported in part by the US NSF, ARO, ONR, and the Sloan Foundation. AI assistance (ChatGPT 5.6) was used in the preparation of this work, including generating code and figures, suggesting and checking proofs, and revising the article. The authors wrote the exposition and prose. All proofs were verified by the authors and we take responsibility for any errors.
\bibliography{main}