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.
110,370 characters
Double machine learning and design in batch adaptive experiments
\maketitle
\begin{abstract}
We consider an experiment with at least two stages or batches and $O(N)$ subjects per batch.
First, we propose a semiparametric treatment effect estimator that efficiently pools information across the batches,
and show it asymptotically dominates alternatives that aggregate single batch estimates.
Then,
we consider the design problem of learning propensity scores for assigning treatment in the later batches of the experiment
to maximize the asymptotic precision of this estimator.
For two common causal estimands,
we estimate this precision using observations from previous batches,
and then solve a finite-dimensional concave maximization problem
to adaptively learn flexible propensity scores that converge to suitably defined optima in each batch at rate $O_p(N^{-1/4})$.
By extending the framework of double machine learning,
we show this rate suffices for our pooled estimator
to attain the targeted precision after each batch,
as long as nuisance function estimates converge at rate $o_p(N^{-1/4})$.
These relatively weak rate requirements enable the investigator to avoid the common practice of discretizing the covariate space for design and estimation in batch adaptive experiments
while maintaining the advantages of pooling.
Our numerical study shows that
such discretization often leads to substantial asymptotic and finite sample precision losses
outweighing any gains from design.
\end{abstract}
\section{Introduction}\label{sec:introduction}
In sequential experimentation, we can use
earlier observations to adjust
our treatment allocation policy for subsequent observations and thereby gain improved
estimation of causal effects in the overall study.
For instance, for an experiment with one treatment arm and one control arm,~\citet{neyman1934} showed that choosing the number of subjects in each arm to be proportional to the outcome standard deviation of that arm
minimizes the variance of the treatment effect estimate based on the difference in means.
While these standard deviations are unknown,
they can be estimated using the initial data.
Then the Neyman allocation can be approximated to improve the sample efficiency of the remainder of the experiment~\citep{hahn2011adaptive,blackwell2022batch, zhao2023adaptive,dai2023clip}.
We study a version of this design problem
for an experiment
divided into a small number of stages or \emph{batches}.
The design, or treatment assignment mechanism,
can be updated \emph{adaptively} for later batches based on the observations from earlier batches
to improve
the precision of the causal estimate computed at the end of the experiment.
A salient feature of our setting is knowledge of pre-treatment covariates that can further improve precision.
Thus, we conceptualize our design problem as choosing a \emph{propensity score} for each batch.
The propensity score specifies the probability that a subject receives treatment
given their covariates
(throughout, we consider the setting of a binary treatment).
The propensity score is well known to be a key mathematical object to be estimated in the causal analysis of observational data (e.g.~\citet{rosenbaum1983central}).
In a randomized experiment it is known and under the control of the investigator.
Hence, it can be exploited for design.
Our specific design objective is to minimize an appropriate scalarization of the asymptotic covariance matrix of an estimator that efficiently pools information across all batches of the experiment.
After describing our mathematical notation and setup in Section~\ref{sec:setup},
we present and study an oracle version of this pooled estimator in Section~\ref{sec:pooled_estimation}.
That oracle is given knowledge of some nuisance parameters,
typically infinite-dimensional mean or variance functions.
It pertains to
a so-called
``non-adaptive batch experiment"
where treatment is assigned with possibly varying but nonrandom propensity scores across batches.
In certain cases,
our oracle pooled estimator asymptotically dominates the best possible
alternative that aggregates single batch estimates,
regardless of the per-batch propensities.
This justifies designing for the pooled estimator instead of an
aggregation-based alternative.
For such a design procedure to be useful,
however, we must show that the targeted asymptotic precision of the oracle pooled estimator is in fact attainable in a batched experiment where the propensity scores used are adaptive (data-dependent)
and the nuisance parameters need to be estimated.
We address these challenges in Section~\ref{sec:batch_clt}
by extending the framework of double machine learning formalized by~\citet{chernozhukov2018double},
hereafter DML,
to what we call a ``convergent split batch adaptive experiment," or CSBAE.
In a CSBAE,
observations in each batch are split into $K$ folds,
and treatment is assigned in each fold according to an adaptive propensity score that only depends on observations from previous batches within the same fold.
Within a given batch,
the $K$ adaptive propensity scores are further required to converge to a common limit
at rate $O_p(N^{-1/4})$ in root mean square (RMS).
Our DML extension then shows that by plugging in estimated nuisance functions,
we can construct a feasible estimator in a CSBAE
that is asymptotically equivalent to
the oracle pooled estimator computed on the limiting non-adaptive batch experiment.
The nuisance function estimates only need to converge at the rate $o_p(N^{-1/4})$.
Section~\ref{sec:batch_learning} details a finite dimensional concave maximization procedure (Algorithm~\ref{alg:csbae})
that provably constructs a CSBAE
for which the limiting propensities in each batch are sequentially optimal within a function class satisfying standard complexity conditions.
Hence, we can effectively design for our pooled estimator in a batch adaptive experiment
with theoretical guarantees that our final estimator will indeed attain the targeted optimal asymptotic precision,
even with a fairly flexible propensity score learning method and nonparametric machine learning estimates for nuisance functions.
To the best of our knowledge,
existing work either designs for a less efficient alternative to our pooled estimator,
or discretizes the covariate space at the design and estimation stages to construct a feasible variant of the pooled estimator on a batch adaptive experiment.
Our simulations in Section~\ref{sec:simulations} suggest that the latter approach in particular
can lead to substantial precision losses that swamp any gains from design,
even with a moderate number of continuous covariates.
Thus, for the practitioner,
we provide an end-to-end design and estimation procedure to efficiently handle continuous covariates in a batched experiment.
\subsection{Related work}
\label{sec:lit_review}
There has been substantial research interest in adaptive experiment designs in recent years.
In many applications,
treatment assignments are updated in an attempt to maximize the (expected)
response values
of either those in the experiment,
as in adaptive bandit algorithms~\citep{russo2018tutorial,hao2020adaptive},
or those in the superpopulation from which the experimental subjects are assumed to arrive~\citep{xu2018fully,kasy2021adaptive}.
Inference on data collected from these algorithms can be challenging since the treatment assignment rules often do not converge~\citep{hadad2021confidence,zhang2020inference,zhang2021statistical}.
By contrast,
in our setting where the goal is purely statistical (maximizing asymptotic precision of the treatment estimate),
the design objective is a static propensity score to be learned consistently.
An interesting direction for further study would be to design for a mixture of both statistical and non-statistical objectives.
For example, one might expand the literature on tie-breaker designs~\citep{owen2020optimizing,morrison2022optimality,li2023general,kluger2023kernel}) to the setting of batched experiments.
The present work can be viewed as an extension of~\citet{hahn2011adaptive} in several directions.
Those authors considered a two batch experiment
to estimate the average treatment effect (ATE) as precisely as possible.
Using data from the first batch to estimate variance functions,
they estimate the asymptotic variance of a pooled version of the semiparametric efficient ATE estimator of~\citet{hirano2003efficient}
for a coarsely discretized covariate.
Then,
they learn a propensity score for the second batch that approximately minimizes this variance.
The covariate discretization ensures
nuisance functions and optimal propensities can be estimated at parametric $O_p(N^{-1/2})$ rates without parametric assumptions.
Consequently,
a feasible version of the pooled estimator
indeed attains the targeted asymptotic variance on the batch adaptive experiment.
We generalize this pooling construction beyond the setting of ATE estimation and relax these rate requirements
to those described in the previous section.
This permits more efficient handling of continuous covariates through nonparametric nuisance function estimates
and more flexible adaptive propensity scores.
Other approaches to extend the work of~\citet{hahn2011adaptive} include~\citet{kato2020efficient},
who consider an online setting where subjects from a stationary superpopulation enter one at a time,
without batches.
Similarly, the literature on covariate-adjusted response-adaptive (CARA) designs has focused on different but related objectives, both statistical and ethical~\citep{zhang2007asymptotic,zhu2023covariate}.
In the batched setting,
~\citet{tabord-meehan2022stratification}
proposes a method to learn a variance-minimizing stratification of the covariate space of fixed size,
avoiding the need to discretize the space prior to observing the data as in~\citet{hahn2011adaptive}.~\citet{cytrynbaum2021designing} showed that by performing a form of highly stratified treatment assignment called local randomization,
consistent variance function estimates from the first batch
make it possible to attain the semiparametric lower bound for ATE estimation in the second batch with optimal propensity score
without having to estimate the conditional mean functions.
Neither of these approaches, however,
maintains the efficiency advantages of pooling.
They also do not immediately extend beyond ATE estimation.
\section{Setup and notation}
\label{sec:setup}
Let $T \geqslant 2$ be the number of batches in the experiment.
Each subject $i=1,\ldots,N_t$ in batch $t=1,\ldots,T$
has observed covariates $X_{ti}\in\mathbb{R}^d$ and potential
outcomes $Y_{ti}(0),Y_{ti}(1)\in\mathbb{R}$.
We place these in the vectors $S_{ti}=(X_{ti}^\top, Y_{ti}(0),Y_{ti}(1))^\top$ which
are exogenous in our model.
Let $Z_{ti} \in \{0,1\}$ be the binary treatment indicator for this subject,
which is controlled by the investigator.
Under the usual stable unit value treatment assumption (SUTVA),
the observed outcome is
\begin{equation}
\label{eq:sutva}
Y_{ti} = Z_{ti}Y_{ti}(1) + (1-Z_{ti})Y_{ti}(0).
\end{equation}
Then the available data for the subject is $W_{ti}=(X_{ti}^\top,Z_{ti},Y_{ti})^\top \in \mathcal{W}$.
We assume that the vectors $S_{ti}\stackrel{\mathrm{iid}}{\sim} P^S$ for $t=1,\ldots,T$ and $i=1,\ldots,N_t$,
for some distribution $P^S$.
Appendix~\ref{app:nonstationary} relaxes this assumption to
permit certain forms of non-stationarity across batches,
such as covariate shifts.
It will be convenient to define the functions
\begin{equation}
\label{eq:cond_moments}
m_0(z,x)=\mathbb{E}[Y(z) \!\mid\! X=x]
\quad \text{and}\quad
v_0(z,x)=\textnormal{Var}(Y(z) \!\mid\! X=x),
\end{equation}
for $z\in\{0,1\}$ and $x\in\mathcal{X}$.
These expectations are taken under $P^S$.
Let $N = N_1 + \ldots +N_T$.
Then as in~\citet{hahn2011adaptive},~\citet{che2023adaptive} and others,
we consider a proportional asymptotic regime
\begin{equation}
\label{eq:prop_asymp_limit}
\lim_{N \rightarrow \infty} \frac{N_t}{N} = \kappa_t \in (0,1), \quad t=1,\ldots,T
\end{equation}
as $N \rightarrow \infty$.
In settings where the batch sizes are fully controlled by the experimenter,
it may be theoretically preferable to make initial batch sample sizes a vanishing fraction of the total sample size~\citep{zhao2023adaptive}.
However,
in many settings
the batch sizes are exogenously constrained to satisfy~\eqref{eq:prop_asymp_limit}
unless observations are discarded.
For various $q \geqslant 1$ and probability measures $P$ on some space $\Omega$,
it will be useful to consider function norms of the form
\[
\|f\|_{q,P} = \left(\int |f(w)|^q \mathrm{d} P(w)\right)^{1/q}
\]
for $f\in L^q(P)$.
We will
use propensity scores
denoted by $e(\cdot)$ with various subscripts.
A propensity score $e(\cdot)$ specifies $e(x)=\textnormal{Pr}(Z=1 \mid X=x)$,
the probability of treatment conditional on covariates.
We will typically require propensity scores to lie in $\mathcal{F}_{\gamma}$,
the set of all measurable functions on $\mathcal{X}$ taking on values in the interval $[\gamma,1-\gamma]$ for some $\gamma \in [0,1/2)$.
We use $\|A\|$ to denote the square root of the sum of the squared entries of any vector, matrix, or tensor $A$.
For any integer $p \geqslant 1$,
$\mathbb{S}_+^p$ will denote the set of symmetric positive semidefinite $p \times p$ real matrices,
and $\mathbb{S}_{++}^p$ will be the set of symmetric positive definite $p \times p$ real matrices.
Finally, for any real vector $v$, we write $v^{\otimes 2}=vv^\top$.
We summarize the preceding requirements for the data generating process in Assumption~\ref{assump:DGP}.
Assumption~\ref{assump:DGP} does not impose any restrictions on the treatment assignment process,
which will be discussed at length in subsequent sections.
\begin{assumption}[Data generating process]
\label{assump:DGP}
For some fixed number of batches $T \geqslant 2$,
the vectors
\[
S_{ti}=(X_{ti},Y_{ti}(0),Y_{ti}(1)), \quad 1 \leqslant t \leqslant T,\quad 1 \leqslant i \leqslant N_t
\]
are independent and identically distributed (i.i.d.) from a distribution $P^S$.
Furthermore, the sample sizes $N_t$ satisfy~\eqref{eq:prop_asymp_limit},
and the vector $W_{ti}=(X_{ti},Z_{ti},Y_{ti})$ is observed
where the outcomes $Y_{ti}$ satisfy the SUTVA assumption~\eqref{eq:sutva}.
\end{assumption}
\subsection{Estimands and score equations}
Consider the setting where $T=1$ (so we can drop the batch subscript $t$),
and the observations $W_1,\ldots,W_N$
are i.i.d.
Suppose additionally that~\eqref{eq:sutva} holds
along with the unconfoundedness assumption
\[
(Y_i(0),Y_i(1)) \perp \!\!\! \perp Z_i \mid X_i, \quad i=1,\ldots,N.
\]
Then many popular causal estimands $\theta_0 \in \Theta \subseteq \mathbb{R}^p$ are identified by a score equation
\[
\mathbb{E}[s(W;\theta_0,\nu_0,e_0)]=0.
\]
In this score equation,
$\nu_0$ is a vector of possibly infinite-dimensional nuisance parameters lying in a nuisance set $\mathcal{N}$,
and $e_0=e_0(\cdot):\mathcal{X} \rightarrow [0,1]$ is the propensity score.
Following Section 3.1 of~\citet{chernozhukov2018double},
we will assume for simplicity that the score $s(\cdot)$ is \emph{linear} in the sense that
\begin{equation}
\label{eq:linear_score}
s(w;\theta,\nu,e) = s_a(w;\nu,e)\theta + s_b(w;\nu,e), \quad \forall w \in \mathcal{W},\ \theta \in \Theta,\ (\nu,e) \in \mathcal{N} \times \mathcal{F}_{\gamma}
\end{equation}
for some $\gamma \in [0,1/2)$,
$s_a(\cdot,\nu,e):\mathcal{W} \rightarrow \mathbb{R}^{p \times p}$,
and $s_b(\cdot,\nu,e): \mathcal{W} \rightarrow \mathbb{R}^p$.
When $T>1$,
propensity scores may vary across batches by design or external constraints.
For any propensity $e=e(\cdot)$ and integrable function $f:\mathcal{W} \rightarrow \mathbb{R}$,
we use the subscripted notation $\mathbb{E}_e[f(W)]=\int f(w) \mathrm{d} P_e(w)$ where $P_e=P_e^W$ is the distribution of
$W=(X,Z,Y)=(X,Z,ZY(1)+(1-Z)Y(0))$
induced by $S=(X,Y(0),Y(1)) \sim P^S$
and $Z \!\mid\! X \sim \textnormal{Bern}(e(X))$ under the SUTVA assumption~\eqref{eq:sutva}.
Further let $P^X$ be the marginal distribution of $X$ under $S \sim P^S$.
Then we will require the following score equations to hold for some $\gamma \in [0,1/2)$
to identify the our causal estimand $\theta_0$:
\begin{equation}
\label{eq:score}
\mathbb{E}_e[s(W;\theta_0,\nu_0,e')]=0, \quad \forall e,e' \in \mathcal{F}_{\gamma}.
\end{equation}
Note that~\eqref{eq:score}
requires the score $s(\cdot)$ to have mean 0 when \emph{any} propensity score $e'(\cdot) \in \mathcal{F}_{\gamma}$ is plugged in.
This plug-in propensity $e'(\cdot)$ may differ from the propensity $e(\cdot) \in \mathcal{F}_{\gamma}$ used for treatment assignment
in the experiment that generated the observations $W$.
Such a robustness property is satisfied by definition so long as the score $s(\cdot)$ is \emph{doubly robust}.
It is required to ensure the validity of the pooled estimator that we propose in Section~\ref{sec:pooled_estimation}.
That estimator requires plugging in a mixture propensity score that averages propensities across all batches $t=1,\ldots,T$.
While identification of $\theta_0$ within each batch is possible by only requiring~\eqref{eq:score} to hold when $e(\cdot)=e'(\cdot)$,
this will not be sufficient to ensure validity of our pooled estimator.
We formally restate our requirements on identification of the estimand $\theta_0$ in Assumption~\ref{assump:identification}.
\begin{assumption}[Estimand identification]
\label{assump:identification}
The estimand $\theta_0 \in \mathbb{R}^p$ of interest satisfies~\eqref{eq:score} for some $\gamma \in [0,1/2)$,
some nuisance parameters $\nu_0$ lying in a known convex set $\mathcal{N}$,
and some score $s(\cdot)$ satisfying~\eqref{eq:linear_score}.
\end{assumption}
The first estimand we are motivated by is the ATE, given by
$\theta_{0,\textnormal{ATE}}=\mathbb{E}[Y(1)-Y(0)]\in\mathbb{R}$
in our notation.
For an investigator interested in modeling how the treatment effect varies with $X$,
they may instead wish to estimate the regression parameter $\theta_{0,\textnormal{PL}} \in \mathbb{R}^p$ under a linear treatment effect assumption
\begin{equation}
\label{eq:pl_assumption}
\mathbb{E}[Y(1)-Y(0) \mid X] = \psi(X)^{\top}\theta_{0,\textnormal{PL}}.
\end{equation}
See~\citet{robinson1988root} for background on semiparametric estimation of $\theta_{0,\textnormal{PL}}$ under~\eqref{eq:pl_assumption},
which characterizes the well-known ``partially linear model."
We show next that both $\theta_{0,\textnormal{ATE}}$ and $\theta_{0,\textnormal{PL}}$ are identified by score functions
that are linear in the sense of~\eqref{eq:linear_score}
and robust in the sense of~\eqref{eq:score},
and hence identifiable according to Assumption~\ref{assump:identification}.
\begin{example}[ATE estimation]
\label{ex:aipw_score}
Let $\theta_0=\theta_{0,\textnormal{ATE}}$ be the estimand of interest.
Now consider the augmented inverse propensity weighting (AIPW) score function
\begin{equation}
\label{eq:aipw_score}
s_{\textnormal{AIPW}}(W;\theta,\nu,e) = m(1,X)-m(0,X)+\frac{Z(Y-m(1,X))}{e(X)}-\frac{(1-Z)(Y-m(0,X))}{1-e(X)}-\theta
\end{equation}
for nuisance parameter $\nu=(m(0,\cdot),m(1,\cdot))$.
For each $\gamma > 0$,
it is well known that $\mathbb{E}_{e}[s_{\textnormal{AIPW}}(W;\theta_0,\nu_0,e')]=0$ for any $e(\cdot),e'(\cdot)$ in $\mathcal{F}_{\gamma}$
when $\nu_0=\nu_{0,\textnormal{AIPW}}=(m_0(0,\cdot), m_0(1,\cdot))$ lies in
the nuisance set $\mathcal{N}=\mathcal{N}_{\textnormal{AIPW}}=L^1(P^X) \times L^1(P^X)$.
Hence,
$s_{\textnormal{AIPW}}(\cdot)$ satisfies the score equation~\eqref{eq:score}.
This score is also linear
because $s_{\textnormal{AIPW}}=s_{\textnormal{AIPW},a}\theta+s_{\textnormal{AIPW},b}$ for
\begin{align*}
s_{\textnormal{AIPW},a}(W;\nu,e) & = -1,\quad\text{and} \\
s_{\textnormal{AIPW},b}(W;\nu,e) & = m(1,X)-m(0,X)+\frac{Z(Y-m(1,X))}{e(X)}-\frac{(1-Z)(Y-m(0,X))}{1-e(X)}.
\end{align*}
which completes the task of showing that $\theta_{0,\textnormal{ATE}}$ satisfies Assumption~\ref{assump:identification}.
\end{example}
\begin{example}[Partially linear model]
\label{ex:epl_score}
Suppose the linear treatment effect assumption~\eqref{eq:pl_assumption} holds and $\theta_0=\theta_{0,\textnormal{PL}}$ is the estimand of interest. Now consider the weighted least squares score
\begin{equation}
\label{eq:epl_score}
s_{\textnormal{EPL}}(W;\theta,\nu,e) = w(X;\nu,e)(Z-e(X))(Y-m(0,X)-Z\psi(X)^{\top}\theta)\psi(X)
\end{equation}
for nonnegative weights
$$w(X,\nu,e)= (v(0,x)e(x)+v(1,x)(1-e(x)))^{-1}$$
with nuisance parameter $\nu=(m(0,\cdot),v(0,\cdot),v(1,\cdot)^\top$.
Let $\mathcal{F}(\mathcal{X};I)$ be the set of all measurable functions $f:\mathcal{X} \rightarrow I$.
Then if $\nu_0=(m_0(0,\cdot),v_0(0,\cdot),v_0(1,\cdot))$ lies in the nuisance set $\mathcal{N}=\mathcal{N}_{\textnormal{EPL}}=L^2(P^X) \times \mathcal{F}(\mathcal{X};[c,\infty)) \times \mathcal{F}(\mathcal{X};[c,\infty))$
for some $c>0$,
we have $\mathbb{E}_e[s_{\textnormal{EPL}}(W;\theta_0,\nu_0,e')]=0$ for any $e(\cdot),e'(\cdot)$ in $\mathcal{F}_0$.
Furthermore,
$s_{\textnormal{EPL}}(\cdot)$ is linear,
because $s_{\textnormal{EPL}}=s_{\textnormal{EPL},a}\theta+s_{\textnormal{EPL},b}$ for
\begin{align*}
s_{\textnormal{EPL},a}(W;\nu,e) & = -w(X;\nu,e)Z(Z-e(X))\psi(X)\psi(X)^{\top},\quad\text{and} \\
s_{\textnormal{EPL},b}(W;\nu,e) &= w(X;\nu,e)(Z-e(X))(Y-m(0,X))\psi(X).
\end{align*}
Thus, $\theta_{0,\textnormal{EPL}}$ satisfies Assumption~\ref{assump:identification} with $\gamma=0$ and the score $s_{\textnormal{EPL}}(\cdot)$.
Note that the score equations~\eqref{eq:score} hold for any nonnegative weight functions $w(\cdot,\nu,e) \in L^1(P^X)$,
though the specific choice in $s_{\textnormal{EPL}}(\cdot)$ is semiparametrically efficient~\citep{chamberlain1992efficiency,ma2006efficient}.
\end{example}
Some of our later results pertain to general estimands identified by Assumption~\ref{assump:identification}.
Others will be specialized to the settings of Examples~\ref{ex:aipw_score} and~\ref{ex:epl_score}.
\section{Oracle pooled and aggregate estimation (non-adaptive)}
\label{sec:pooled_estimation}
Here we propose and analyze an oracle estimator $\hat{\theta}^*$ of a generic estimand $\theta_0$ satisfying Assumption~\ref{assump:identification} with score $s(\cdot)$.
It is an oracle in the sense that it uses the unknown
true value of the nuisance parameter $\nu_0$ from the score equations~\eqref{eq:score}.
We discuss feasible estimation of $\theta_0$, including estimation of $\nu_0$,
in Section~\ref{sec:batch_clt}.
The estimator $\hat{\theta}^*$ pools observations across all batches $t=1,\ldots,T$.
We then prove a central limit theorem (CLT) for it
in the setting of a \emph{non-adaptive batch experiment}
where treatment in each batch $t=1,\ldots,T$ is assigned according to a fixed (non-random) propensity score $e_t(\cdot)$:
\begin{definition}[Non-adaptive batch experiment]
\label{def:non_adaptive_batch_experiment}
A \textbf{non-adaptive batch experiment} has data generating process satisfying Assumption~\ref{assump:DGP} and treatment assignments satisfying
\[
Z_{ti} = \bm{1}(U_{ti} \leqslant e_t(X_{ti})), \quad t=1,\ldots,T,\ i=1,\ldots,N_t
\]
for some \emph{nonrandom} (i.e., non-adaptive) batch propensity scores $e_1(\cdot),\ldots,e_T(\cdot)$ and uniformly distributed random variables $\{U_{ti} \mid t=1,\ldots,T,i=1,\ldots,N_t\}$ that are i.i.d.\ and independent of the vectors $\{S_{ti} \mid t=1,\ldots,T,i=1,\ldots,N_t\}$.
\end{definition}
Next, we compare the pooled estimator $\hat{\theta}^*$ to an alternative oracle
that makes an optimal linear aggregation of per-batch
estimates.
It also satisfies a CLT,
but we show that our pooling
strategy dominates aggregation in terms of efficiency in the setting of Examples~\ref{ex:aipw_score} and~\ref{ex:epl_score}.
We remark that some authors
(e.g.~\citet{tabord-meehan2022stratification}) refer to this aggregation approach as pooling,
but we reserve that term for pooling data, not estimators.
\subsection{Pooled oracle estimator}
The main idea behind the construction of our oracle estimator $\hat{\theta}^*$
is as follows.
After collecting the observations from all batches $t=1,\ldots,T$ in a non-adaptive batch experiment,
we ignore the batch structure and pool together the observations across batches.
Now consider a random draw $W=(X,Z,Y)$ from these pooled observations $\{W_{ti} \mid 1 \leqslant t \leqslant T, 1 \leqslant i \leqslant N_t\}$.
In the notation of Section~\ref{sec:setup},
it is straightforward to show that the distribution of $W$ is $P_{e_{0,N}}=\sum_{t=1}^T (N_t/N)P_{e_t}$,
where $e_{0,N}(\cdot)$ is the mixture propensity score
\begin{equation}
\label{eq:e_0_N}
e_{0,N}(x) = \textnormal{Pr}(Z=1 \mid X) = \sum_{t=1}^T \frac{N_t}{N}e_t(x).
\end{equation}
It will also be helpful to define the limiting mixture propensity score $e_0(\cdot)$ under the proportional asymptotics~\eqref{eq:prop_asymp_limit}:
\begin{equation}
\label{eq:e_0}
e_0(x) = \sum_{t=1}^T \kappa_te_t(x).
\end{equation}
When a particular set of nonrandom batch propensities $e_1(\cdot),\ldots,e_T(\cdot)$ is relevant,
we omit an additional subscript by letting
$\mathbb{E}_{0,N}[f(W)]=\mathbb{E}_{e_{0,N}}[f(W)]$
and $\mathbb{E}_0[f(W)]=\mathbb{E}_{e_0}[f(W)]$.
Using this notation, the oracle pooled estimator $\hat{\theta}^*$ is derived by solving the sample analogue of the mixture score equations $\mathbb{E}_{0,N}[s(W;\theta_0,\nu_0,e_{0,N})]=0$ for $\theta$:
\begin{equation}
\label{eq:theta_star}
\hat{\theta}^* = -\biggl(\frac{1}{N}\sum_{t=1}^T\sum_{i=1}^{N_t} s_a(W_{ti};\nu_0,e_{0,N})\biggr)^{-1} \biggl(\frac{1}{N}\sum_{t=1}^T\sum_{i=1}^{N_t} s_b(W_{ti};\nu_0,e_{0,N})\biggr).
\end{equation}
This $\hat\theta^*$ is only defined when the matrix inverse
in~\eqref{eq:theta_star} exists.
That will be the case with probability tending to 1 so long as $\mathbb{E}_0[s_a(W;\nu_0,e_0)]$ is invertible,
which we will require in our theoretical results.
One such result is a CLT for $\hat{\theta}^*$.
The proofs of all our technical results are provided in Appendix~\ref{app:proofs}.
\begin{proposition}[Oracle CLT]
\label{prop:oracle_clt}
Let $\theta_0 \in \mathbb{R}^p$ be an estimand satisfying Assumption~\ref{assump:identification} for some $\gamma \in [0,1/2)$,
some score $s(\cdot)$,
and some nuisance parameters $\nu_0$.
Suppose observations $\{W_{ti} \mid t=1,\ldots,T,i=1,\ldots,N_t\}$ are collected from a non-adaptive batch experiment with batch propensities $e_1(\cdot),\ldots,e_T(\cdot) \in \mathcal{F}_{\gamma}$,
and define $e_{0,N}=e_{0,N}(\cdot)$ and $e_0=e_0(\cdot)$ as in~\eqref{eq:e_0_N} and~\eqref{eq:e_0}, respectively.
Further assume the following conditions hold:
\begin{enumerate}
\item \label{cond:score_continuity} For some sequence $\delta_N \downarrow 0$, we have
\begin{align*}
\bigl(\mathbb{E}_0[\|s_a(W;\nu_0,e_{0,N})-s_a(W;\nu_0,e_0)\|^2]\bigr)^{1/2} & \leqslant \delta_N,\quad\text{and} \\
\bigl(\mathbb{E}_0[\|s(W;\theta_0,\nu_0,e_{0,N})-s(W;\theta_0,\nu_0,e_0)\|^2]\bigr)^{1/2} & \leqslant \delta_N.
\end{align*}
\item \label{cond:s_a_invertibility} $\mathbb{E}_0[s_a(W;\nu_0,e_0)]$ is invertible and $\mathbb{E}_0[\|s(W;\theta_0,\nu_0,e_0)\|^2] < \infty$.
\item \label{cond:moment_boundedness} For some $q>2$ and $C < \infty$ we have $\mathbb{E}_0[\|s(W;\theta_0,\nu_0,e_{0,N})\|^q]\leqslant C$ for all sufficiently large $N$.
\end{enumerate}
Then with $\hat{\theta}^*$ as defined in~\eqref{eq:theta_star}, we have
\begin{align*}
\sqrt{N}(\hat{\theta}^*-\theta_0) & \stackrel{\mathrm{d\,\,}}\rightarrow \mathcal{N}(0, V_0)
\end{align*}
where
\begin{align*}
V_0 & = \bigl(\mathbb{E}_0[s_a(W;\nu_0,e_0)]\bigr)^{-1}\bigl(\mathbb{E}_0[s(W;\theta_0,\nu_0,e_0)^{\otimes 2}]\bigr)\bigl(\mathbb{E}_0[s_a(W;\nu_0,e_0)]\bigr)^{-1}.
\end{align*}
\end{proposition}
\begin{proof}
See Appendix~\ref{proof:prop:oracle_clt}.
\end{proof}
For the estimands
$\theta_{0,\textnormal{ATE}}$ and $\theta_{0,\textnormal{PL}}$,
following~\eqref{eq:theta_star}
we use the scores $s_{\textnormal{AIPW}}(\cdot)$ and $s_{\textnormal{EPL}}(\cdot)$,
respectively to derive the oracle estimates
\begin{align}
\hat{\theta}^*_{\textnormal{AIPW}} & = \frac{1}{N} \sum_{t=1}^T\sum_{i=1}^{N_t} s_{\textnormal{AIPW},b}(W_{ti};\nu_{0,\textnormal{AIPW}},e_{0,N}) \quad \text{ (recall $s_{\textnormal{AIPW}, a}=-1$)},\quad \text{and} \label{eq:theta_star_aipw} \\
\hat{\theta}^*_{\textnormal{EPL}} & = -\biggl(\frac{1}{N} \sum_{t=1}^T\sum_{i=1}^{N_t} s_{\textnormal{EPL},a}(W_{ti};\nu_{0,\textnormal{EPL}},e_{0,N})\biggr)^{-1}\biggl(\frac{1}{N}\sum_{t=1}^T\sum_{i=1}^{N_t} s_{\textnormal{EPL},b}(W_{ti};\nu_{0,\textnormal{EPL}},e_{0,N})\biggr), \label{eq:theta_star_epl}
\end{align}
We now specialize the generic oracle CLT of Proposition~\ref{prop:oracle_clt} to these two estimators under some regularity conditions.
\begin{assumption}
\label{assump:ate_regularity}[Regularity for estimating $\theta_{0,\textnormal{ATE}}$]
For some $C< \infty$ and $q>2$,
we have $(\mathbb{E}[|Y(z)|^q])^{1/q} \leqslant C$ and $\mathbb{E}[Y(z)^2 \mid X=x] \leqslant C$ for all $z=0,1$ and $x \in \mathcal{X}$.
\end{assumption}
\begin{assumption}
\label{assump:pl_regularity}[Regularity for estimating $\theta_{0,\textnormal{PL}}$]
For some $C< \infty$ and $q>2$,
Assumption~\ref{assump:ate_regularity} holds.
Additionally, $\|\psi(x)\| \leqslant C$ for all $x \in \mathcal{X}$,
and there exists $c>0$ such that $v_0(z,x) \geqslant c$ for all $z=0,1$ and $x \in \mathcal{X}$.
Finally, the linear treatment effect assumption~\eqref{eq:pl_assumption} holds.
\end{assumption}
\begin{corollary}[Oracle CLT for $\hat{\theta}_{\textnormal{AIPW}}^*$]
\label{cor:ate_oracle_clt}
Suppose Assumption~\ref{assump:ate_regularity} holds,
and let $\{W_{ti} \mid t=1,\ldots,T,i=1,\ldots,N_t\}$ be observations from a non-adaptive batch experiment with batch propensities $e_1(\cdot),\ldots,e_T(\cdot) \in \mathcal{F}_{\gamma}$ for some $\gamma>0$.
Then $\sqrt{N}(\hat{\theta}^*_{\textnormal{AIPW}}-\theta_{0,\textnormal{ATE}}) \stackrel{\mathrm{d\,\,}}\rightarrow \mathcal{N}(0,V_{0,\textnormal{AIPW}})$
where
\begin{align}
V_{0,\textnormal{AIPW}} & = \mathbb{E}\left[\frac{v_0(1,X)}{e_0(X)} + \frac{v_0(0,X)}{1-e_0(X)} + (m_0(1,X)-m_0(0,X)-\theta_{0,\textnormal{ATE}})^2\right]. \label{eq:v_0_aipw}
\end{align}
\end{corollary}
\begin{proof}
See Appendix \ref{proof:cor:ate_oracle_clt}.
\end{proof}
\begin{corollary}[Oracle CLT for $\hat{\theta}_{\textnormal{EPL}}^*$]
\label{cor:pl_oracle_clt}
Suppose Assumption~\ref{assump:pl_regularity} holds,
and let $\{W_{ti} \mid t=1,\ldots,T,i=1,\ldots,N_t\}$ be observations from a non-adaptive batch experiment with batch propensities $e_1(\cdot),\ldots,e_T(\cdot) \in \mathcal{F}_0$,
where $\mathbb{E}[e_0^2(X)(1-e_0(X))^2\psi(X)\psi(X)^{\top}] \in \mathbb{S}^p_{++}$.
Then $\sqrt{N}(\hat{\theta}^*_{\textnormal{EPL}}-\theta_{0,\textnormal{PL}}) \stackrel{\mathrm{d\,\,}}\rightarrow \mathcal{N}(0,V_{0,\textnormal{EPL}})$,
where
\begin{align}
V_{0,\textnormal{EPL}} & = \left(\mathbb{E}\left[\frac{e_0(X)(1-e_0(X))}{v_0(0,X)e_0(X) + v_0(1,X)(1-e_0(X))}\psi(X)\psi(X)^{\top}\right]\right)^{-1}. \label{eq:v_0_pl}
\end{align}
\end{corollary}
\begin{proof}
See Appendix \ref{proof:cor:pl_oracle_clt}.
\end{proof}
\subsection{The aggregated oracle estimator}
When the covariate space $\mathcal{X}$ is finite,
the oracle pooled estimator $\hat{\theta}^*_{\textnormal{AIPW}}$ is equivalent to an oracle variant of the estimator proposed by~\citet{hahn2011adaptive}.
The complexity of constructing a \emph{feasible} pooled estimator
for an \emph{adaptive} experiment in more general settings
has led other authors to instead consider single batch estimators that can lose considerable efficiency.
For instance,~\citet{cytrynbaum2021designing} proposes simply discarding the first batch in a two-batch experiment when computing the final estimate,
which is clearly inadmissible in our proportional asymptotic regime~\eqref{eq:prop_asymp_conv}.
Section 3.2 of~\citet{tabord-meehan2022stratification} suggests instead taking a linear aggregation of estimates computed separately on each batch,
as described above.
We will now show that even the best linearly aggregated oracle estimator is asymptotically dominated by our pooled estimators $\hat{\theta}^*_{\textnormal{AIPW}}$ and $\hat{\theta}^*_{\textnormal{EPL}}$.
While~\citet{tabord-meehan2022stratification} hypothesized in their Appendix C.2 that this may be true for ATE estimation,
they do not pursue this further as are unable to construct a feasible pooled estimator
attaining the targeted oracle variance
when batch propensities are chosen adaptively using their stratification trees.
By contrast, our design approach
will allow us to construct such an estimator using our extension of double machine learning in Section~\ref{sec:batch_clt}.
Fix a non-adaptive batch experiment with batch propensities $e_1(\cdot),\ldots,e_T(\cdot)$.
For each batch $t=1,\ldots,T$,
consider an (oracle) estimator $\hat{\theta}_{t,\textnormal{AIPW}}^*$ for $\theta_{0,\textnormal{ATE}}$ computed by solving the empirical analogue of the score equations $\mathbb{E}_{e_t}[s_{\textnormal{AIPW}}(W;\theta_{0,\textnormal{ATE}},\nu_{0,\textnormal{AIPW}},e_t)]=0$
that averages only those observations in batch $t$:
\[
\hat{\theta}_{t,\textnormal{AIPW}}^* = -\biggl(\frac{1}{N_t}\sum_{i=1}^{N_t} s_{\textnormal{AIPW},a}(W_{ti};\nu_{0,\textnormal{AIPW}},e_t)\biggr)^{-1} \biggl(\frac{1}{N_t}\sum_{i=1}^{N_t} s_{\textnormal{AIPW},b}(W_{ti};\nu_{0,\textnormal{AIPW}},e_t)\biggr).
\]
By applying Proposition~\ref{prop:oracle_clt} with a single batch,
for each $t=1,\ldots,T$ we obtain the CLT
\begin{align*}
\sqrt{N_t}(\hat{\theta}_{t,\textnormal{AIPW}}^*-\theta_{0,\textnormal{ATE}}) & \stackrel{\mathrm{d\,\,}}\rightarrow \mathcal{N}(0,V_{t,\textnormal{AIPW}}),
\end{align*}
where
\begin{align*}
V_{t,\textnormal{AIPW}} &=A_{t,\textnormal{AIPW}}^{-1}B_{t,\textnormal{AIPW}}A_{t,\textnormal{AIPW}}^{-1},\quad\text{for} \\
A_{t,\textnormal{AIPW}} & = \mathbb{E}_{e_t}[s_{\textnormal{AIPW},a}(W;\nu_{0,\textnormal{AIPW}},e_t)],\quad\text{and} \\
B_{t,\textnormal{AIPW}} & = \mathbb{E}_{e_t}[s(W;\theta_{0,\textnormal{ATE}},\nu_{0,\textnormal{AIPW}},e_t)^{\otimes 2}].
\end{align*}
Now, as stated in~\citet{slud2018combining},
the asymptotically unbiased linear combination of $\hat{\theta}_{1,\textnormal{AIPW}}^*,\ldots,\hat{\theta}_{T,\textnormal{AIPW}}^*$ with the smallest asymptotic covariance matrix with respect to the semidefinite ordering is the inverse covariance weighted estimator
\[
\hat{\theta}_{\textnormal{AIPW}}^{*,(\mathrm{LA})} = \biggl(\,\sum_{t=1}^T \kappa_t V_{t,\textnormal{AIPW}}^{-1}\biggr)^{-1} \sum_{t=1}^T \kappa_t V_{t,\textnormal{AIPW}}^{-1}\hat{\theta}_{t,\textnormal{AIPW}}^*.
\]
This optimal linearly aggregated estimator $\hat{\theta}_{\textnormal{AIPW}}^{*,(\mathrm{LA})}$ satisfies the CLT
\[
\sqrt{N}\bigl(\hat{\theta}_{\textnormal{AIPW}}^{*,(\mathrm{LA})}-\theta_{0,\textnormal{ATE}}\bigr) \stackrel{\mathrm{d\,\,}}\rightarrow \mathcal{N}(0,V_{\textnormal{AIPW}}^{(\mathrm{LA})}), \quad V_{\textnormal{AIPW}}^{(\mathrm{LA})} = \biggl(\sum_{t=1}^T \kappa_t V_{t,\textnormal{AIPW}}^{-1}\biggr)^{-1}.
\]
We can similarly define the linearly aggregated oracle estimator $\hat{\theta}_{\textnormal{EPL}}^{*,(\mathrm{LA})}$ for $\theta_{0,\textnormal{PL}}$ based on combining per-batch estimates from the score $s_{\textnormal{EPL}}(\cdot)$,
which satisfies a CLT $\sqrt{N}(\hat{\theta}^{*,(\mathrm{LA})}_{\textnormal{EPL}}-\theta_{0,\textnormal{PL}}) \stackrel{\mathrm{d\,\,}}\rightarrow \mathcal{N}(0,V_{\textnormal{EPL}}^{(\mathrm{LA})})$.
Here $V_{\textnormal{EPL}}^{(\mathrm{LA})}$ is given by replacing
$s_{\textnormal{AIPW}}(\cdot)$ with $s_{\textnormal{EPL}}(\cdot)$, $\nu_{0,\textnormal{AIPW}}$ with $\nu_{0,\textnormal{EPL}}$,
and $\theta_{0,\textnormal{ATE}}$ with $\theta_{0,\textnormal{PL}}$ in the definition of $V_{\textnormal{AIPW}}^{(\mathrm{LA})}$.
Our main result is then that regardless of the batch propensities,
$\hat{\theta}_{\textnormal{AIPW}}^{*,(\mathrm{LA})}$ and $\hat{\theta}_{\textnormal{EPL}}^{*,(\mathrm{LA})}$ are asymptotically dominated by our pooled estimators $\hat{\theta}_{0,\textnormal{AIPW}}$ and $\hat{\theta}_{0,\textnormal{EPL}}$, respectively.
This motivates our work in Section~\ref{sec:batch_learning} that designs for these estimators.
\begin{theorem}[Pooling dominates linear aggregation]
\label{thm:pooled_covariance}
Under the conditions of Corollary~\ref{cor:ate_oracle_clt},
$V_{0,\textnormal{AIPW}} \leqslant V^{(\mathrm{LA})}_{\textnormal{AIPW}}$.
Under the conditions of Corollary~\ref{cor:pl_oracle_clt},
$V_{0,\textnormal{EPL}} \preccurlyeq V^{(\mathrm{LA})}_{\textnormal{EPL}}$.
\end{theorem}
\begin{proof}
See Appendix~\ref{proof:thm:pooled_covariance}.
\end{proof}
\section{Feasible pooled estimation in batch adaptive experiments}
\label{sec:batch_clt}
The oracle estimator $\hat\theta^*$ of~\eqref{eq:theta_star} depends on nuisance parameters $\nu_0$ that are unknown in practice.
Additionally, recall that our CLT for $\hat\theta^*$ (Proposition~\ref{prop:oracle_clt}) holds for a non-adaptive batch experiment.
Our goal is to choose propensities adaptively in each batch to improve precision.
Therefore, we
would like to
develop a feasible estimator $\hat\theta$ that attains the targeted asymptotic variance for experiments where treatment is assigned \emph{adaptively},
even when the nuisance parameters $\nu_0$ must be estimated.
As mentioned above,
our construction of such a feasible estimator $\hat{\theta}$ is based on extending the double machine learning (DML) framework of~\citet{chernozhukov2018double}.
The main requirements for $\hat{\theta}$ to have the same asymptotic variance as the corresponding oracle are convergence rate guarantees for both nuisance parameter estimates and the adaptive propensities.
Our DML extension ensures that these rate requirements can be made sub-parametric,
enabling the use of somewhat flexible machine learning methods.
The typical DML setting assumes access to a single sample $W_1,\ldots,W_N$ of i.i.d.\ observations.
An example of this setting is a non-adaptive batch experiment with $T=1$ and propensity $e_0(\cdot)$.
Then a standard DML estimator is based on two ingredients: a Neyman orthogonal score and cross-fitting.
Neyman orthogonality of the score $s(\cdot)$ at $(\nu_0,e_0)$ means a local insensitivity of the score equations to perturbations in $(\nu_0,e_0)$ in any direction:
\[
\frac{\partial}{\partial \lambda} \mathbb{E}_{e_0}\bigl[s(W_i;\theta_0,\nu_0+\lambda(\nu-\nu_0),e_0+\lambda(e-e_0))\bigr]=0, \quad \forall (\nu,e) \in \mathcal{N} \times \mathcal{F}_{\gamma}.
\]
It is well known that the scores $s_{\textnormal{AIPW}}(\cdot)$ and $s_{\textnormal{EPL}}(\cdot)$ are Neyman orthogonal~\citep{chernozhukov2018double}.
Given a Neyman orthogonal score $s(\cdot)$,
DML proceeds by constructing an estimator $\hat{\theta}$ by cross-fitting.
In cross-fitting,
the indices $1,\ldots,N$ are partitioned into $K$ (roughly) equally sized folds $\mathcal{I}_1,\ldots,\mathcal{I}_K$.
Then $\hat{\theta}$ is computed as the solution to the empirical score equations
$N^{-1} \sum_{k=1}^K \sum_{i \in \mathcal{I}_k} s(W_i;\hat{\theta},\hat{\nu}^{(-k)},\hat{e}^{(-k)}) = 0$,
where for each $k=1,\ldots,K$,
$\hat{\nu}^{(-k)}$ and $\hat{e}^{(-k)}(\cdot)$ are estimates of $\nu_0$ and $e_0(\cdot)$,
respectively.
Each pair $(\hat{\nu}^{(-k)},\hat{e}^{(-k)}(\cdot))$ depends only on the observations $\{W_i \mid i \notin \mathcal{I}_k\}$ outside fold $k$.
The sample splitting ensures that for each $k=1,\ldots,K$,
the estimates $(\hat{\nu}^{(-k)},\hat{e}^{(-k)})$ are independent of the observations in fold $k$.
By the arguments of~\citet{chernozhukov2018double},
such independence is key to guarantee that the feasible estimator $\hat{\theta}$ is equivalent (up to first order asymptotics) to the oracle $\hat{\theta}^*$ solving
$N^{-1} \sum_{k=1}^K \sum_{i \in \mathcal{I}_k} s(W_i;\hat{\theta},\nu_0,e_0) = 0$ (cf.~\eqref{eq:theta_star}),
even when the estimates $(\hat{\nu}^{(-k)},\hat{e}^{(-k)})$ converge at sub-parametric rates.
To maintain this independence in an \emph{adaptive} batched experiment,
we require sample splitting at the \emph{design} stage,
as illustrated in Figure~\ref{fig:csbae},
along with convergence of the adaptive propensity scores.
Our notion of a convergent split batch adaptive experiment (CSBAE) in Definition~\ref{def:CSBAE} formalizes this.
The main idea is to split the observations in \emph{every} batch $t=1,\ldots,T$ into $K$ folds.
Re-using the notation above from the standard DML setting with $T=1$,
we let $\mathcal{I}_k$ denote the set of batch and observation indices $(t,i)$ assigned to fold $k=1,\ldots,K$.
Then the adaptive propensity used to assign treatment to a subject in batches $t=2,\ldots,T$
is allowed to only depend on observations in previous batches from the same fold as this subject.
To ensure that the adaptivity does not introduce any additional variability
(up to first-order asymptotics)
into the final estimator,
a CSBAE requires these adaptive propensities to converge to nonrandom limits $e_1(\cdot),\ldots,e_T(\cdot)$ at RMS rate $O_p(N^{-1/4})$.
While this convergence requirement may appear restrictive,
in Section~\ref{sec:batch_learning} we show how it can be ensured by design by solving an appropriate finite-dimensional concave maximization procedure.
Moreover, the limiting propensity scores from this procedure will be provably optimal,
in a sense we make more precise in Section~\ref{sec:batch_learning}.
\begin{figure}[!ht]
\includegraphics[width=\linewidth]{SBAE.pdf}
\caption{A graphical representation of the dependencies among the observations in two batches of a CSBAE (Definition~\ref{def:CSBAE}).
Note the lack of vertical arrows,
indicating independence between the observations in fold $j$ and fold $k$ for any distinct $j$, $k$ in $\{1,\ldots,K\}$.}
\label{fig:csbae}
\end{figure}
\begin{definition}[Convergent split batch adaptive experiment]
\label{def:CSBAE}
A \textbf{convergent split batch adaptive experiment (CSBAE)}
is an experiment with data generating process satisfying Assumption~\ref{assump:DGP}
where each observation index $(t,i)$ is assigned
to one of $K$ folds $\mathcal{I}_1,\ldots,\mathcal{I}_K$.
The fold assignments are such that $n_{t,k}= |\{(t,i) \in \mathcal{I}_k \mid i=1,\ldots,N_t\}|$,
the number of observations in batch $t$ assigned to fold $k$,
satisfies $|n_{t,k}-N_t/K| \leqslant 1$ for all $t=1,\ldots,T$, $k=1,\ldots,K$.
Now let $P_{N,t}^{X,(k)}$ be the empirical distribution on $\{X_{ti} \mid (t,i) \in \mathcal{I}_k, 1 \leqslant i \leqslant N_t\}$,
the covariates in batch $t$ and fold $k$,
and define $\mathcal{S}_t^{X,(k)}$
to be the $\sigma$-algebra generated by the covariates $\{X_{ti} \mid (t,i) \in \mathcal{I}_k, 1 \leqslant i \leqslant N_t\}$ in batch $t$ and fold $k$ along with the observations $\{W_{ui} \mid (u,i) \in \mathcal{I}_k, u=1,\ldots,t-1,i=1,\ldots,N_u\}$ in fold $k$ and any of the previous batches $1,\ldots,t-1$.
We further require the following for each batch $t=1,\ldots,T$ and fold $k=1,\ldots,K$:
\begin{enumerate}
\item Treatment is assigned according to an adaptive propensity $\hat{e}_t^{(k)}(\cdot)$ that is measurable with respect to $\mathcal{S}_t^{X,(k)}$.
That is, the treatment indicators can be represented as
\begin{equation}
\label{eq:CSBAE_treatment}
Z_{ti}=\bm{1}(U_{ti} \leqslant \hat{e}_t^{(k)}(X_{it})), \quad (t,i) \in \mathcal{I}_k
\end{equation}
where $\{U_{ti}:1 \leqslant t \leqslant T, 1 \leqslant i \leqslant N_t\}$ is a collection of i.i.d.\ uniformly distributed random variables independent of the vectors $\{S_{ti} \mid t=1,\ldots,T,i=1,\ldots,N_t\}$.
\item For some nonrandom propensity $e_t(\cdot)$,
the adaptive propensity $\hat{e}_t^{(k)}(\cdot)$ satisfies
\begin{equation}
\label{eq:csbae_limit}
\|\hat{e}_t^{(k)}-e_t\|_{2,P_{N,t}^{X,(k)}} = O_p(N^{-1/4}), \quad t=1,\ldots,T, \ k=1,\ldots,K.
\end{equation}
\end{enumerate}
\begin{remark}
The left-hand side of equation~\eqref{eq:csbae_limit} uses an $L^2$ norm on the empirical distribution $P_{N,t}^{X,(k)}$ of the covariates of the subjects that will be assigned treatment according to the learned propensity $\hat{e}_t^{(k)}(\cdot)$.
These covariates will also be used to learn $\hat{e}_t^{(k)}(\cdot)$ itself in our propensity learning procedure of Section~\ref{sec:batch_learning}.
Thus, we can interpret~\eqref{eq:csbae_limit} as a rate requirement on the ``in-sample" convergence of $\hat e_t^{(k)}(\cdot)$.
\end{remark}
\end{definition}
Given a CSBAE,
our
feasible estimator is
\begin{equation}
\label{eq:theta_hat}
\hat{\theta} = -\biggl(\frac{1}{N} \sum_{k=1}^K \sum_{(t,i) \in \mathcal{I}_k} s_a(W_{ti};\hat{\nu}^{(-k)},\hat{e}^{(-k)})\biggr)^{-1}\biggl(\frac{1}{N} \sum_{k=1}^K \sum_{(t,i) \in \mathcal{I}_k} s_b(W_{ti};\hat{\nu}^{(-k)},\hat{e}^{(-k)})\biggr).
\end{equation}
As in the standard (single batch) DML setting,
for each $k=1,\ldots,K$,
$\hat{\nu}^{(-k)}$ and $\hat{e}^{(-k)}(\cdot)$ are estimates of the nuisance parameters $\nu_0$ and the mixture propensity $e_{0,N}(\cdot)$ defined in~\eqref{eq:e_0_N}, respectively,
that depend only on the observations $\{W_{ti} \mid (t,i) \notin \mathcal{I}_k\}$ outside fold $k$.
These observations are fully independent of the observations in fold $k$ (across \emph{all} batches $t=1,\ldots,T$)
by the construction of a CSBAE.
As in the single batch case,
given $o_p(N^{-1/4})$ convergence of the estimators $\hat{\nu}^{(-k)}$ to $\nu_0$,
this independence along with~\eqref{eq:csbae_limit} enable a DML-style argument that $\hat{\theta}$ is asymptotically equivalent to the oracle $\hat{\theta}^*$ under~\eqref{eq:csbae_limit} computed on a counterfactual
non-adaptive batch experiment with propensities $e_1(\cdot),\ldots,e_T(\cdot)$.
This argument proceeds by coupling the treatment indicators $Z_{ti}$ in the
CSBAE with counterfactual treatment indicators $\tilde{Z}_{ti}=\bm{1}(U_{ti} \leqslant e_t(X_{ti}))$.
\begin{assumption}[Score properties and convergence rates for estimating nuisance parameters and the mixture propensity in a CSBAE]
\label{assump:dml}
Observations $\{W_{ti}:1 \leqslant t \leqslant T, 1 \leqslant i \leqslant N_t\}$ are collected from a CSBAE with limiting batch propensities $e_1(\cdot),\ldots,e_T(\cdot)$.
Additionally,
the estimand $\theta_0$ of interest is identified as in Assumption~\ref{assump:identification} by some score $s(\cdot)$,
nuisance parameters $\nu_0 \in \mathcal{N}$,
and $\gamma \in [0,1/2)$,
such that the propensity collection $\mathcal{F}_{\gamma}$ contains $e_1(\cdot),\ldots,e_T(\cdot)$.
Defining $W_{ti}(z)=(Y_{ti}(z),X_{ti},z)$ for $z=0,1$,
the score $s(\cdot)$ has the following properties:
\begin{enumerate}[label=(\alph*)]
\item \label{cond:regularity} The matrix $\mathbb{E}_0[s_a(W;\nu_0,e_0)]$ is invertible and $\mathbb{E}_0\bigl[\|s(W;\theta_0,\nu_0,e_0)\|^2\bigr] < \infty$.
\item \label{cond:score_differentiability} The mapping $\lambda \mapsto \mathbb{E}_{0,N}[s(W;\theta_0,\nu_0+\lambda(\nu-\nu_0),e_{0,N}+\lambda(e-e_{0,N})]$ is twice continuously differentiable on $[0,1]$ for each $(\nu,e) \in \mathcal{T}=\mathcal{N} \times \mathcal{F}_{\gamma}$.
\item \label{cond:potential_score_equations}
All propensities $e(\cdot) \in \mathcal{F}_{\gamma}$ satisfy
\begin{equation}
\label{eq:potential_score_diff}
\mathbb{E}\bigl[s(W_{ti}(1);\theta_0,\nu_0,e)-s(W_{ti}(0);\theta_0,\nu_0,e) \!\mid\! X_{ti}\bigr] = 0,
\end{equation}
for all $t=1,\dots,T$ and $i=1,\dots,N_t$.
\end{enumerate}
Also, there exist estimators $\hat{\nu}^{(-k)}$ and $\hat{e}^{(-k)}(\cdot)$ of $\nu_0$ and the mixture propensity $e_{0,N}(\cdot)$ defined in~\eqref{eq:e_0_N},
respectively,
that depend only on the observations outside fold $k$ of the CSBAE.
Next,
there are nonrandom subsets $\mathcal{T}_N \subseteq \mathcal{T}$ containing $(\nu_0,e_{0,N}(\cdot))$,
such that for all $k=1,\ldots,K$,
$\textnormal{Pr}((\hat{\nu}^{(-k)},\hat{e}^{(-k)}(\cdot)) \in \mathcal{T}_N) \rightarrow 1$ as $N \rightarrow \infty$.
The sets $\mathcal{T}_N$ shrink quickly enough for the following to hold for all $(\nu,e)\in\mathcal{T}_N$,
all $\lambda\in(0,1)$
and all $z\in\{0,1\}$
when $N$ is sufficiently large:
\begin{align}
\bigg\|\frac{\partial}{\partial \lambda} \mathbb{E}_{0,N}[s(W;\theta_0,\nu_0+\lambda(\nu-\nu_0),e_{0,N}+\lambda(e-e_{0,N}))] \Big|_{\lambda=0} \bigg\| & \leqslant N^{-1/2}\delta_N \label{eq:neyman_orthogonality} \\
\bigg\|\frac{\partial^2}{\partial \lambda^2} \mathbb{E}_{0,N}[s(W;\theta_0,\nu_0+\lambda(\nu-\nu_0),e_{0,N}+\lambda(e-e_{0,N}))] \bigg\| & \leqslant N^{-1/2}\delta_N
\label{eq:vanishing_second_derivatives} \\
\left(\mathbb{E}_{0}\bigl[\|s_a(W;\nu,e)-s_a(W;\nu_0,e_0)\|^2\bigr]\right)^{1/2} & \leqslant \delta_N \label{eq:s_a_consistency} \\
\left(\mathbb{E}_{0}\bigl[\|s(W;\theta_0,\nu,e)-s(W;\theta_0,\nu_0,e_0)\|^2\bigr]\right)^{1/2} & \leqslant \delta_N
\label{eq:s_consistency} \\
\bigl(\mathbb{E}_{0}\bigl[\|s_a(W(z);\nu,e)\|^q\bigr]\bigr)^{1/q} & \leqslant C
\label{eq:s_a_moment_condition} \\
\bigl(\mathbb{E}_{0}\bigr[\|s(W(z);\theta_0,\nu,e\|^q\bigr]\bigr)^{1/q} & \leqslant C.
\label{eq:s_moment_condition}
\end{align}
Finally,
letting $\mathcal{S}^{(-k)}$ be the $\sigma$-algebra generated by the observations $\{W_{ti}:(t,i) \notin \mathcal{I}_k\}$ outside fold $k$ across all batches 1 through $T$,
we require
\begin{equation}
\label{eq:S_t_neg_k}
S_t^{(-k)}(z) = o_p(N^{-1/4}), \quad z \in \{0,1\}, \quad k=1,\ldots,K
\end{equation}
where
\[
S_t^{(-k)}(z) = \sqrt{
\frac{1}{n_{t,k}} \sum_{i:(t,i) \in \mathcal{I}_k} \Bigl\|\mathbb{E}\bigl[s(W_{ti}(z);\hat{\nu}^{(-k)},\hat{e}^{(-k)})-s(W_{ti}(z);\nu_0,e_{0,N}) \bigm| \mathcal{S}^{(-k)},X_{ti}\bigr]\Bigr\|^2}.
\]
\end{assumption}
\begin{theorem}[Feasible CLT for a CSBAE]
\label{thm:batch_clt}
Suppose Assumption~\ref{assump:dml} holds. Then for $\hat{\theta}$ defined in~\eqref{eq:theta_hat}
there exists a non-adaptive batch experiment with propensities $e_1(\cdot),\ldots,e_T(\cdot)$
for which
\[
\hat{\theta} = \hat{\theta}^* + o_p(N^{-1/2}).
\]
Here $\hat{\theta}^*$ is the oracle~\eqref{eq:theta_star} computed on this non-adaptive batch experiment.
Then $\sqrt{N}(\hat{\theta}-\theta_0) \stackrel{\mathrm{d\,\,}}\rightarrow \mathcal{N}(0,V_0)$ where $V_0$ is the limiting covariance derived in Proposition~\ref{prop:oracle_clt}.
\end{theorem}
\begin{proof} See Appendix~\ref{proof:thm:batch_clt}.
\end{proof}
In Assumption~\ref{assump:dml},
the conditions~\ref{cond:regularity} and~\ref{cond:score_differentiability}
along with the inequalities~\eqref{eq:neyman_orthogonality} through~\eqref{eq:s_moment_condition} are direct extensions of Assumptions 3.1 and 3.2 in~\citet{chernozhukov2018double}
for ordinary DML ($T=1$).
The equations~\eqref{eq:potential_score_diff} and~\eqref{eq:S_t_neg_k} are additional requirements that enable the dependence across batches in a CSBAE to be sufficiently weak
so that $\hat{\theta}$ computed on the CSBAE is asymptotically equivalent to a version of it computed on the limiting non-adaptive batch experiment.
We can show that Assumption~\ref{assump:dml} is satisfied for estimating $\theta_{0,\textnormal{ATE}}$ and $\theta_{0,\textnormal{PL}}$ with $s_{\textnormal{AIPW}}(\cdot)$ and $s_{\textnormal{EPL}}(\cdot)$, respectively
under simple rate conditions on nuisance parameter and propensity estimation rates that mirror those in Section 5 of~\citet{chernozhukov2018double} for single-batch DML estimators.
Then we apply Theorem~\ref{thm:batch_clt} to construct feasible pooled estimators $\hat{\theta}_{\textnormal{AIPW}}$ and $\hat{\theta}_{\textnormal{EPL}}$ as special cases of equation~\eqref{eq:theta_hat},
which attain the oracle asymptotic variances $V_{0,\textnormal{AIPW}}$ and $V_{0,\textnormal{EPL}}$ defined in~\eqref{eq:v_0_aipw} and~\eqref{eq:v_0_pl}, respectively.
\begin{corollary}[Feasible pooled estimation of $\theta_{0,\textnormal{ATE}}$ in a CSBAE]
\label{cor:AIPW}
Suppose observations are collected from a CSBAE
for which the regularity conditions of Assumption~\ref{assump:ate_regularity} hold
for some $q>2$ and $C<\infty$
and the limiting batch propensities $e_1(\cdot),\ldots,e_T(\cdot)$ are in $\mathcal{F}_{\gamma}$
for some $\gamma > 0$.
Additionally, for each $k=1,\ldots,K$,
suppose we have estimates $\hat{m}^{(-k)}(\cdot)$ and $\hat{e}^{(-k)}(\cdot)$ of the mean function $m_0(\cdot)$ and mixture propensity $e_{0,N}(\cdot)$,
respectively, both depending only on the observations outside fold $k$,
such that the following are true:
\begin{enumerate}
\item $\|\hat{m}^{(-k)}(z,\cdot)-m_0(z,\cdot)\|_{2,P^X}=o_p(N^{-1/4}), \quad z=0,1,$
\item $\|\hat{e}^{(-k)}-e_{0,N}\|_{2,P^X} = O_p(N^{-1/4}),$
\item $\|\hat{m}^{(-k)}(z,\cdot)-m_0(z,\cdot)\|_{q,P^X} \leqslant C, \quad z=0,1$,
\item $\hat{e}^{(-k)}(\cdot) \in \mathcal{F}_{\gamma}$ with probability tending to 1.
\end{enumerate}
Then $N^{1/2}(\hat{\theta}_{\textnormal{AIPW}}-\theta_{0,\textnormal{ATE}}) \stackrel{\mathrm{d\,\,}}\rightarrow \mathcal{N}(0,V_{0,\textnormal{AIPW}})$.
\end{corollary}
\begin{proof}
See Appendix \ref{proof:cor:AIPW}.
\end{proof}
\begin{corollary}[Feasible estimation of $\theta_{0,\textnormal{PL}}$ in a CSBAE]
\label{cor:pl}
Suppose observations are collected from a CSBAE
for which the regularity conditions of Assumption~\ref{assump:pl_regularity} hold
for some $q>2$ and $0<c<C<\infty$,
and for which the limiting batch propensities $e_1(\cdot),\ldots,e_T(\cdot)$ are in $\mathcal{F}_0$
with $\mathbb{E}[e_0^2(X)(1-e_0(X))^2\psi(X)\psi(X)^{\top}] \in \mathbb{S}^p_{++}$.
For each $k=1,\ldots,K$,
assume we have estimates $\hat{m}^{(-k)}(0,\cdot)$, $\hat{v}^{(-k)}(\cdot,\cdot)$, and $\hat{e}^{(-k)}(\cdot)$ of the mean function $m_0(0,\cdot)$,
the variance function $v_0(\cdot, \cdot)$,
and the mixture propensity $e_{0,N}(\cdot)$,
respectively,
all depending only on the observations outside fold $k$,
such that the following are true:
\begin{enumerate}
\item $\|\hat{m}^{(-k)}(0,\cdot)-m_0(0,\cdot)\|_{2,P^X}=o_p(N^{-1/4})$,
\item $\|\hat{v}^{(-k)}(z,\cdot)-v_0(z,\cdot)\|_{2,P^X}=o_p(1), \quad z=0,1$,
\item $\|\hat{e}^{(-k)}-e_{0,N}\|_{2,P^X} = O_p(N^{-1/4})$,
\item $\|\hat{m}^{(-k)}(0,\cdot)-m_0(0,\cdot)\|_{q,P^X} \leqslant C,$
\item $\inf_{x \in \mathcal{X}} \hat{v}^{(-k)}(z,x) \geqslant c, \quad z=0,1$.
\end{enumerate}
Then $N^{1/2}(\hat{\theta}_{\textnormal{EPL}}-\theta_{0,\textnormal{PL}}) \stackrel{\mathrm{d\,\,}}\rightarrow \mathcal{N}(0,V_{0,\textnormal{EPL}})$.
\end{corollary}
\begin{proof}
See Appendix \ref{proof:cor:pl}.
\end{proof}
We now compare the rate requirements in Corollaries~\ref{cor:AIPW} and~\ref{cor:pl}
with those needed to prove a feasible CLT for the linearly aggregated estimators discussed in Section~\ref{sec:pooled_estimation}.
Consider a batch adaptive experiment \emph{without} sample splitting,
so that we assign treatment using the propensities $\hat{e}_1(\cdot),\ldots,\hat{e}_T(\cdot)$
where each $\hat{e}_t(\cdot)$ is possibly random but can only depend on the observations in batches $1,\ldots,t-1$,
and converges to some nonrandom $e_t(\cdot)$
at \emph{any} rate;
in particular, this rate may be slower than the $O_p(N^{-1/4})$ rate of~\eqref{eq:csbae_limit}.
Then the standard DML results in Section 5 of~\citet{chernozhukov2018double} show that we can construct feasible estimators $\hat{\theta}_{t,\textnormal{AIPW}}$ and $\hat{\theta}_{t,\textnormal{EPL}}$ that are asymptotically equivalent to the oracle single-batch estimators $\hat{\theta}_{t,\textnormal{AIPW}}^*$ and $\hat{\theta}_{t,\textnormal{EPL}}^*$ in Section~\ref{sec:pooled_estimation},
so long as we plug in the true propensity score $\hat{e}_t(\cdot)$ used to assign treatment
and use cross-fitting at the \emph{estimation} stage,
even if all rates in Corollaries~\ref{cor:AIPW} and~\ref{cor:pl} are weakened to $o_p(1)$.
Unfortunately,
we cannot extend this construction to our pooled estimator $\hat{\theta}$ when $T \geqslant 2$,
since $\hat{\theta}$ plugs in estimates of the mixture propensity $e_{0,N}(\cdot)$,
which does not correspond (in general) to a propensity actually used for treatment in any batch.
Our numerical studies in Section~\ref{sec:simulations} suggest,
however, that the weaker rate requirements for a feasible CLT for the linearly aggregated estimators
(compared to our pooled estimators) make little difference in practice.
Indeed,
there we also see finite sample advantages to pooling,
beyond those predicted by the asymptotics.
\section{Batch adaptive learning of the optimal propensity score}
\label{sec:batch_learning}
We now discuss how to learn adaptive propensity scores $\hat{e}_t^{(k)}(\cdot)$ that satisfy~\eqref{eq:csbae_limit}
with limiting propensity scores $e_t(\cdot)$ that maximize asymptotic precision of the final estimators $\hat{\theta}_{\textnormal{AIPW}}$ and $\hat{\theta}_{\textnormal{EPL}}$ constructed in the previous section,
as measured by their asymptotic covariance matrices $V_{0,\textnormal{AIPW}}$ and $V_{0,\textnormal{EPL}}$.
This generates a CSBAE on which,
by Theorem~\ref{thm:batch_clt},
feasible estimators $\hat{\theta}_{0,\textnormal{AIPW}}$ and $\hat{\theta}_{0,\textnormal{EPL}}$ achieve the targeted asymptotic variances $V_{0,\textnormal{AIPW}}$ and $V_{0,\textnormal{EPL}}$.
While $V_{0,\textnormal{AIPW}}$ is scalar
and so there is no ambiguity in what it means to maximize asymptotic precision of $\hat{\theta}_{\textnormal{AIPW}}$,
when $\theta_{0,\textnormal{PL}}$ is multivariate ($p>1$),
$V_{0,\textnormal{EPL}}$ is a matrix.
To handle the multivariate setting,
we follow classical literature on experiment design in regression models
by scalarizing an \emph{information matrix} $\mathcal{I} \in \mathbb{S}_+^p$ (typically an inverse covariance matrix) using an \textit{information function} $\Psi:\mathbb{S}_+^p \rightarrow \mathbb{R} \cup \{-\infty\}$.
See textbooks such as~\citet{atkinson2007optimum,pukelsheim2006optimal} for more background.
We generically write $\mathcal{I}=\mathcal{I}(e,\eta)$ to emphasize that in our setting,
$\mathcal{I}$ will be indexed by a propensity score $e=e(\cdot)$ in some function class $\mathcal{F}_*$
along with some unknown nuisance functions $\eta=\eta(\cdot)$ to be estimated.
Then the design objective is to learn an optimal propensity $e^*(\cdot)$,
in the sense of maximizing $\Psi(\mathcal{I})$:
\begin{equation}
\label{eq:e_star}
e^*(\cdot) \in \operatorname*{arg\,max}_{e \in \mathcal{F}_*} \Psi(\mathcal{I}(e;\eta)).
\end{equation}
The nuisance functions $\eta$ in the information matrix are distinct from the nuisance functions $\nu$ in the score function.
Examples of $\eta$ for information matrices based on $V_{0,\textnormal{AIPW}}$ and $V_{0,\textnormal{EPL}}$ are given in Section~\ref{sec:apply_concave_maximization}.
\subsection{Generic convergence rates with concave maximization}
We start by considering a feasible procedure to learn the propensity $e^*(\cdot)$ of~\eqref{eq:e_star} when the information matrix $\mathcal{I}$
and the function class $\mathcal{F}_*$ being optimized over generically satisfy Assumption~\ref{assump:info_matrix} below.
An explanation of how to apply this generic procedure to construct an appropriate CSBAE that designs for the estimators $\hat{\theta}_{\textnormal{AIPW}}$ and $\hat{\theta}_{\textnormal{EPL}}$ is deferred to
Section~\ref{sec:apply_concave_maximization}.
Here we instead focus on a precise exposition of our propensity learning procedure and the technical assumptions needed to ensure convergence,
without explicitly invoking any notation or setup from previous sections.
We begin with some generic structure on the information matrix $\mathcal{I}$ and the function class $\mathcal{F}_*$.
Roughly,
we require the information matrix to be strongly concave in $e(\cdot)$
and the function class $\mathcal{F}_*$ to be not too complex.
We also allow for interval constraints on the per-batch expected proportion of subjects treated.
\begin{assumption}[Generic optimization setup]
\label{assump:info_matrix}
The information matrix $\mathcal{I}=\mathcal{I}(e,\eta)$ in~\eqref{eq:e_star} takes the form
\begin{equation}
\label{eq:generic_info}
\mathcal{I}(e;\eta) = \mathbb{E}_{X \sim P}\left[f(e(X),\eta(X))\right]
\end{equation}
where $\eta:\mathcal{X} \rightarrow \mathcal{W}$ is a vector of possibly unknown functions taking values in some compact set $\mathcal{W} \subseteq \mathbb{R}^r$,
and $P$ is a \textbf{covariate distribution} from which i.i.d.\ observations $X_1,\ldots,X_n \in \mathcal{X}$ are drawn.
Additionally,
the function $f:[0,1] \times \mathcal{W} \rightarrow \mathbb{S}_+^p$ appearing in~\eqref{eq:generic_info} satisfies the following properties:
\begin{itemize}
\item There exists an extension of $f(\cdot,\cdot)$ with continuous second partial derivatives to an open neighborhood of $\mathbb{R}^{r+1}$ containing $(\delta,1-\delta) \times \mathcal{W}$ for some $\delta < 0$.
\item For some $c>0$ we have
\begin{equation}
\label{eq:f_strong_concavity}
-f''(e,w) \in \mathbb{S}_+^p, \quad\text{and}\quad \textnormal{tr}(f''(e,w)) \leqslant -c,\quad \forall (e,w) \in [0,1] \times \mathcal{W}
\end{equation}
where $f''(e,w)$ denotes the second partial derivative of $f(\cdot,\cdot)$ with respect to the first argument,
evaluated at $(e,w)$.
\end{itemize}
Next, the collection of propensity scores $\mathcal{F}_*$ to be optimized over
takes the form $\mathcal{F}_* = \mathcal{F}_*(m_L,m_H;P) = \{e \in \mathcal{E} \mid m_L \leqslant \mathbb{E}_P[e(X)] \leqslant m_H\}$,
for some \textbf{base propensity class} $\mathcal{E} \subseteq \mathcal{F}_0$
and some known \textbf{budget constraints} $(m_L,m_H)$ with $0 \leqslant m_L \leqslant m_H \leqslant 1$.
The base propensity class $\mathcal{E}$ is convex and closed in $L^2(P)$
and additionally satisfies the following properties for all $n \geqslant 1$:
\begin{enumerate}
\item \label{cond:finite_entropy} There exists
$C<\infty$ for which
\begin{equation}
\label{eq:finite_entropy}
\int_0^1 \sqrt{\log \mathcal{N}(\epsilon, \mathcal{E}, L^2(P_n))} \,\mathrm{d}\epsilon \leqslant C\quad \forall n\quad \text{w.p.1}
\end{equation}
where $\log \mathcal{N}(\epsilon,\mathcal{E},L^2(P_n))$ is the metric entropy of the base propensity class $\mathcal{E}$ in $L^2(P_n)$,
as defined in Appendix~\ref{app:asymptotics},
and $P_n$ is the empirical distribution on the observations $X_1,\ldots,X_n$.
\item \label{cond:transform_En} Given any $x_1,\ldots,x_n \in \mathcal{X}$,
there exists a convex set $E_n \subseteq [0,1]^n$, possibly depending on $x_1,\ldots,x_n$, such that
\begin{itemize}
\item For every $e \in \mathcal{E}$, we must have $(e(x_1),e(x_2),\ldots,e(x_n)) \in E_n$, and
\item For every $(e_1,e_2,\ldots,e_n) \in E_n$, there exists $e \in \mathcal{E}$ with $e(x_i) = e_i$ for all $i=1,\ldots,n$.
\end{itemize}
\item \label{cond:positivity} There exists $e_*(\cdot) \in \mathcal{F}_*$ such that
$\mathbb{E}_P[f(e_*(X),\eta(X))] \succcurlyeq cI$ for some $c>0$.
\item \label{cond:hilo} There exist $e_L(\cdot),e_H(\cdot) \in \mathcal{E}$ such that $\mathbb{E}_P[e_L(X)] > m_L$ and $\mathbb{E}_P[e_H(X)] < m_H$.
\end{enumerate}
\end{assumption}
\begin{remark}
If~\eqref{eq:f_strong_concavity} only holds for $(e,w)$ in $[\gamma,1-\gamma] \times \mathcal{W}$ for some $\gamma \in (0,1/2)$ (instead of on all of $[0,1] \times \mathcal{W}$) and the base propensity class $\mathcal{E}$ is chosen to lie within $\mathcal{F}_{\gamma}$,
then one can rewrite~\eqref{eq:generic_info} as
\[
\mathcal{I}(e;\eta) = \mathbb{E}_{X \sim P}[\tilde{f}(e(X),\eta(X))], \quad \tilde{f}(e,w) = f(L_{\gamma}(e), w)
\]
where $L_{\gamma}(e)=\gamma+(1-2\gamma)e$ is an invertible linear mapping.
Then the remainder of Assumption~\ref{assump:info_matrix} holds as stated with $f$ replaced by $\tilde{f}$ and the base propensity class $\mathcal{E}$ replaced by $\tilde{\mathcal{E}} = \{L_{\gamma}^{-1}(e(\cdot)) \mid e(\cdot) \in \mathcal{E}\}$,
and so all results below depending on Assumption~\ref{assump:info_matrix} hold by considering optimization over $\tilde{\mathcal{E}}$ instead of $\mathcal{E}$.
\end{remark}
Similar to the idea of empirical risk minimization in supervised learning (ERM, e.g.~\citet{vapnik1991principles} and~\citet{montanari2022universality}),
when Assumption~\ref{assump:info_matrix} is satisfied,
our learning procedure replaces the unknown population expectation appearing in the objective~\eqref{eq:e_star} by a sample average over observations $X_1,\ldots,X_n \stackrel{\mathrm{iid}}{\sim} P$,
where the generic sample size $n$ diverges.
In our designed CSBAE,
these observations will be the covariates of those subjects in a given batch $t=1,\ldots,T$ of a CSBAE within a given fold $k=1,\ldots,K$,
so $n$ will be identified with the quantity $n_{t,k}$ of Definition~\ref{def:CSBAE}.
The main computational procedure is then a finite dimensional optimization over the values of the propensity score at the $n$ points $X_1,\ldots,X_n$:
\begin{equation}
\label{eq:e_hat}
(\hat{e}_1,\ldots,\hat{e}_n) \in \operatorname*{arg\,max}_{(e_1,\ldots,e_n) \in F_n} \Psi\biggl(n^{-1} \sum_{i=1}^n f(e_i,\hat{\eta}(X_i))\biggr).
\end{equation}
Here $\hat{\eta}(\cdot)$ is an estimate of the nuisance function $\eta(\cdot)$ in~\eqref{eq:generic_info},
and the optimization set $F_n=F_n(m_L,m_H) \subseteq \mathbb{R}^n$ is defined by
\begin{equation}
\label{eq:F_n}
F_n(m_L,m_H) = \biggl\{(e_1,\ldots,e_n) \in E_n \mid m_L \leqslant \frac1n \sum_{i=1}^n e_i \leqslant m_H\biggr\} \subseteq \mathbb{R}^n
\end{equation}
where $(m_L,m_H)$ are budget constraints as in Assumption~\ref{assump:info_matrix}
and $E_n$ is as in the numbered condition~\ref{cond:transform_En} of that assumption.
We convert the vector $(\hat{e}_1,\ldots,\hat{e}_n)$ in~\eqref{eq:e_hat} to a propensity score $\hat{e}(\cdot)$
by taking $\hat{e}(\cdot)$ to be any member of the base propensity class $\mathcal{E}$ of Assumption~\ref{assump:info_matrix}
with $\hat{e}(X_i)=\hat{e}_i$ for each $i=1,\ldots,n$.
The existence of such a propensity $\hat{e}(\cdot)$ is guaranteed by the numbered condition~\ref{cond:transform_En} in Assumption~\ref{assump:info_matrix}.
Then as in the ERM literature,
we use empirical process arguments to guarantee the learned propensity $\hat{e}(\cdot)$ converges to $e^*(\cdot)$ at rate $O_p(n^{-1/4})$
under appropriate restrictions on the complexity of the base propensity class $\mathcal{E}$.
Our main such restriction is the finite entropy integral requirement in~\eqref{eq:finite_entropy}.
Examples of function classes satisfying this condition
can be found in the literature on empirical processes.
A sufficient condition is that $\log \mathcal{N}(\epsilon,\mathcal{E},L^2(P_n)) \leqslant K\epsilon^{-2+\delta}$ for some $\delta >0$, some
$K<\infty$,
and all $\epsilon > 0$.
Examples~\ref{example:monotone} through~\ref{example:VC_hull} below show that this requirement is loose enough to admit some relatively rich function classes.
\begin{example}[Monotone]
\label{example:monotone}
Suppose $d=1$ and let $\mathcal{E}$ be the set of nondecreasing functions in $\mathcal{F}_0$.
By Lemma 9.11 of~\citet{kosorok2008introduction},
we know that $\log \mathcal{N}(\epsilon,\mathcal{E},L^2(P_n)) \leqslant K/\epsilon$
for each $\epsilon>0$ and some positive universal constant $K<\infty$.
\end{example}
\begin{example}[Lipschitz]
\label{example:lipschitz}
Again suppose $d=1$ and let $\mathcal{E}$ be the set of $L$-Lipschitz functions in $\mathcal{F}_0$
for some fixed $L>0$.
If $\mathcal{X}$ is a bounded closed interval,
then the discussion preceding Example 5.11 of~\citet{wainwright2019high} shows that $\log \mathcal{N}(\epsilon,\mathcal{E},L^2(P_n)) \leqslant K/\epsilon$
for each $\epsilon>0$ and some positive universal constant
$K<\infty$ (which may depend on $L$).
\end{example}
\begin{example}(VC-subgraph class)
\label{example:VC}
Let $\mathcal{E}$ be any subset of $\mathcal{F}_0$ that is closed and convex in $L^2(P)$,
and whose subgraphs are a Vapnik-Chervonenkis (VC) class,
meaning they have a finite VC dimension $V$.
A special case is a fully parametric class like $\{x \mapsto \theta^{\top}\xi(x) \mid \theta^{\top}\bm{1}_p \leqslant 1, \theta \succcurlyeq 0\}$
where $\xi(x) \in [0,1]^p$ is a known set of basis functions
and $\bm{1}_p := (1,\ldots,1) \in \mathbb{R}^p$.
Then by Theorem 2.6.7 of~\citet{van1996weak},
$\mathcal{N}(\epsilon,\mathcal{E},L^2(P_n)) \leqslant K(1/\epsilon)^{2V-2}$ for some universal $K>0$ depending on the VC dimension $V$ of the subgraphs.
Note that $K$ may depend on $p$.
\end{example}
\begin{example}(Symmetric convex hull of VC-subgraph class)
\label{example:VC_hull}
Let $\mathcal{E}_0$ be a VC-subgraph class of functions.
The symmetric convex hull of $\mathcal{E}_0$ is defined as
\[
\textnormal{sconv}(\mathcal{E}_0) = \biggl\{\,\sum_{i=1}^m \omega_i e_i \Bigm| e_i \in \mathcal{E}_0, \sum_{i=1}^m |\omega_i| \leqslant 1 \biggr\}.
\]
Now suppose $\mathcal{E}$ is contained within $\overline{\textnormal{sconv}(\mathcal{E}_0)}$, the pointwise closure of $\textnormal{sconv}(\mathcal{E}_0)$.
Let $V<\infty$ be the VC dimension of the collection of subgraphs of functions in $\mathcal{E}_0$.
Then by Theorem 2.6.9 of~\citet{van1996weak}
we have $\log \mathcal{N}(\epsilon,\mathcal{E},L^2(P_n)) \leqslant K(1/\epsilon)^{2(1-1/V)}$ for all $\epsilon>0$.
For example, with $\textnormal{expit}(x)=\exp(x)/(1+\exp(x))$ we can take
\begin{equation}
\label{eq:E_conv_hull}
\mathcal{E} = \biggl\{\,\sum_{i=1}^m \omega_i \textnormal{expit}(\theta_i^{\top}\varphi(x)) \Bigm| \omega_i \geqslant 0, \sum_{i=1}^m \omega_i \leqslant 1\biggr\}
\end{equation}
where $\varphi(\cdot)$ is any vector of $p$ real-valued basis functions,
$m$ can be made arbitrarily large,
and $\theta_1,\ldots,\theta_m \in \mathbb{R}^p$ are arbitrary.
This choice of $\mathcal{E}$ is evidently a closed and convex subset of $\textnormal{sconv}(\mathcal{E}_0)$ with $\mathcal{E}_0=\{\textnormal{expit}(\theta^{\top}\varphi(x)) \mid \theta \in \mathbb{R}^p\}$.
Note the collection $\mathcal{E}_0$ is indeed a VC-subgraph class by Lemmas 2.6.15 and 2.6.17 of~\citet{van1996weak},
as each function in $\mathcal{E}_0$ is the composition of the monotone function $\textnormal{expit}(\cdot)$
with the $p$-dimensional vector space of functions $\{x \mapsto \theta^{\top}\varphi(x) \mid \theta \in \mathbb{R}^p\}$.
\end{example}
Next, we construct sets $E_n$ that satisfy condition~\ref{cond:transform_En} of Assumption~\ref{assump:info_matrix}
for the base propensity classes in Examples~\ref{example:monotone} to~\ref{example:VC_hull}.
For the set of monotone functions in one dimension (Example~\ref{example:monotone}) we can take
\[
E_n = \{(e_1,\ldots,e_n) \mid 0 \leqslant e_{\pi(1)} \leqslant e_{\pi(2)} \leqslant \cdots \leqslant e_{\pi(n)} \leqslant 1\}
\]
where $\pi(\cdot)$ is the inverse of the function that maps each $i \in \{1,\ldots,n\}$ to the rank of $x_i$ among $x_1,\ldots,x_n$ (with any ties broken in some deterministic way).
For the set of $L$-Lipschitz functions (Example~\ref{example:lipschitz}) we can take
\[
E_n = \bigl\{(e_1,\ldots,e_n) \bigm| |e_{\pi(i)}-e_{\pi(i-1)}| \leqslant L(x_{\pi(i)}-x_{\pi(i-1)}),\, i=2,\dots,n\bigr\}.
\]
For the parametric class in Example~\ref{example:VC} we can take
\[
E_n = \bigl\{(\theta^{\top}\xi(x_1),\ldots,\theta^{\top}\xi(x_n)) \mid \theta^{\top}\bm{1}_p \leqslant 1, \theta \succcurlyeq 0\bigr\}.
\]
Finally, for the class~\eqref{eq:E_conv_hull} in Example~\ref{example:VC_hull} we take
\[
E_n = \biggl\{\,\sum_{i=1}^m \omega_i\bigl(\textnormal{expit}(\theta_i^{\top}\varphi(x_1)),\ldots,\textnormal{expit}(\theta_i^{\top}\varphi(x_n))\bigr) \Bigm| \omega_i \geqslant 0, \sum_{i=1}^m \omega_i \leqslant 1 \biggr\}.
\]
The convergence rates of $\hat{e}(\cdot)$ to $e^*(\cdot)$ will be proven using strong concavity of the design objective~\eqref{eq:e_star} on the space $\mathcal{F}^*$.
This is ensured by Assumption~\ref{assump:info_matrix} along with the following conditions on the information function $\Psi(\cdot)$:
\begin{assumption}[Information function regularity]
\label{assump:Psi}
The information function $\Psi:\mathbb{S}_+^p \rightarrow \mathbb{R} \cup \{-\infty\}$ is concave, continuous, and nondecreasing with respect to the semidefinite ordering $\succcurlyeq$ on $\mathbb{R}^{p \times p}$ and satisfies the following conditions:
\begin{itemize}
\item[(a)] For every $k>0$, $\inf_{B \succcurlyeq kI} \Psi(B) > \sup_{A \in \mathbb{S}_+^p \setminus \mathbb{S}_{++}^p} \Psi(A) =: \Psi_0$.
\item[(b)] $\Psi(\cdot)$ is
twice continuously differentiable on $\mathbb{S}_{++}^p$,
such that for all $0<k<K$, there exists $C>0$ such that $\|\nabla \Psi(A)-\nabla \Psi(B)\| \leqslant C\|A-B\|$ whenever $KI \succcurlyeq A \succcurlyeq kI$ and $KI \succcurlyeq B \succcurlyeq kI$.
\item[(c)] For every $0 < k < K$,
$kI \preccurlyeq A \preccurlyeq KI$ implies $\tilde{k}I \preccurlyeq \nabla \Psi(A) \preccurlyeq \tilde{K}I$ for some $0 < \tilde{k} < \tilde{K}$.
\item[(d)] For every $K<\infty$ and $\tilde{\Psi}_0 > \Psi_0$, there exists $k > 0$ such that for all $0 \preccurlyeq A \preccurlyeq KI$ with $\Psi(A) \geqslant \tilde{\Psi}_0$, we have $A \succcurlyeq kI$.
\end{itemize}
\end{assumption}
We can show that Assumption~\ref{assump:Psi} is satisfied by two common information functions:
the ``$A$-optimality" function $\Psi_a(\cdot) = -\textnormal{tr}((\cdot)^{-1})$ with $\Psi_a(M) :=-\infty$ whenever $M$ is singular,
and the ``$D$-optimality" function $\Psi_d(\cdot) = \log(\det(\cdot))$.
The $A$-optimality criterion corresponds to minimizing the average (asymptotic) variance of the components of the estimand,
while $D$-optimality corresponds to minimizing the volume of the ellipsoid spanned by the columns of the (asymptotic) covariance matrix.
\begin{lemma}
\label{lemma:Psi_cond}
The information functions $\Psi_d(\cdot)$ and $\Psi_a(\cdot)$ satisfy Assumption~\ref{assump:Psi}.
\end{lemma}
\begin{proof}
See Appendix~\ref{proof:lemma:Psi_cond}.
\end{proof}
There are some common information functions that do not satisfy
Assumption~\ref{assump:Psi}.
For example, the ``$E$-optimality" function $\Psi_e(\cdot) = \lambda_{\min}(\cdot)$,
where $\lambda_{\min}(M)$ refers to the smallest eigenvalue of $M \in \mathbb{R}^{p \times p}$,
is not differentiable.
Similarly the function $\Psi_c(\cdot) = -c^{\top}(\cdot)^{-1}c$ (for some fixed $c \in \mathbb{R}^p$),
corresponding to ``$c$-optimality," does not satisfy condition (c).
We leave open the question of whether the $O_p(n^{-1/4})$
convergence rate of $\hat{e}(\cdot)$ to $e^*(\cdot)$ in Lemma~\ref{lemma:concave_maximization} can be extended to these (and other) information functions using different techniques,
and now prove this rate under Assumptions~\ref{assump:info_matrix} and~\ref{assump:Psi}.
\begin{lemma}[Convergence of generic concave maximization routine]
\label{lemma:concave_maximization}
Suppose Assumption~\ref{assump:info_matrix} holds for some
information matrix $\mathcal{I}$ of the form~\eqref{eq:generic_info},
covariate distribution $P$,
and budget constraints $(m_L,m_H)$.
Further assume that for some sequence $\alpha_n \downarrow 0$,
we have an estimate $\hat{\eta}(\cdot)$ of $\eta(\cdot)$ satisfying
\begin{equation}
\label{eq:eta_n}
\hat{\eta}(x) \in \mathcal{W},\ \forall x \in \mathcal{X}, \quad\text{and}\quad \|\hat{\eta}-\eta\|_{2,P_n} = O_p(\alpha_n)
\end{equation}
where $\eta(\cdot)$ is defined by $\mathcal{I}$ by~\eqref{eq:generic_info}.
Then for any information function $\Psi(\cdot)$ satisfying Assumption~\ref{assump:Psi},
the following statements are true:
\begin{enumerate}
\item There exists an optimal propensity function $e^*(\cdot)$ satisfying~\eqref{eq:e_star} which is unique $P$-almost everywhere.
\item \label{cond:F_n} There exist optimal
finite sample treatment probabilities $(\hat{e}_1,\ldots,\hat{e}_n) \in [0,1]^n$ satisfying~\eqref{eq:e_hat},
where $F_n=F_n(m_L,m_H)$ is defined as in~\eqref{eq:F_n}.
\item For any such optimal probabilities $(\hat{e}_1,\ldots,\hat{e}_n)$,
there exists a propensity score $\hat{e}(\cdot) \in \mathcal{E}$ for which $\hat{e}(X_i) = \hat{e}_i$ for each $i=1,\ldots,n$.
Any such function $\hat{e}(\cdot)$ satisfies both $\|\hat{e}-e^*\|_{2,P}=O_p(n^{-1/4} + \alpha_n)$ and $\|\hat{e}-e^*\|_{2,P_n}=O_p(n^{-1/4}+\alpha_n)$.
\end{enumerate}
\end{lemma}
\begin{proof}
See Appendix \ref{proof:lemma:concave_maximization}.
\end{proof}
\subsection{Convergence of batch adaptive designs }
\label{sec:apply_concave_maximization}
We now leverage Lemma~\ref{lemma:concave_maximization} to develop a procedure
(Algorithm~\ref{alg:csbae})
that can learn adaptive propensities $\hat{e}_t^{(k)}(\cdot)$ with the convergence guarantees~\eqref{eq:csbae_limit} so that when used for treatment assignment,
they lead to a CSBAE
with limiting propensities that optimize objectives of the form~\eqref{eq:e_star}
with information matrices based on $V_{0,\textnormal{AIPW}}^{-1}$ and $V_{0,\textnormal{EPL}}^{-1}$.
This shows we can effectively design for the estimators $\hat{\theta}_{\textnormal{AIPW}}$ and $\hat{\theta}_{\textnormal{EPL}}$.
\begin{algorithm}[tb]
\caption{CSBAE for estimating $\theta_{0,\textnormal{ATE}}$ resp. $\theta_{0,\textnormal{PL}}$}
\label{alg:csbae}
\begin{algorithmic}[1]
\REQUIRE Base propensity class $\mathcal{E}$ satisfying conditions in Assumption~\ref{assump:info_matrix}, initial propensity $e_1(\cdot) \in \mathcal{F}_{\epsilon_1}$ for some $\epsilon_1 > 0$,
information function $\Psi(\cdot)$ satisfying Assumption~\ref{assump:Psi}, number of folds $K \geqslant 2$
\FOR{batch $t=1,\ldots,T$}
\STATE Observe subject covariates $X_{t1},\ldots,X_{tN_t}$
\STATE Split subject indices $i=1,\ldots,N_t$ into $K$ folds $\mathcal{I}_1,\ldots,\mathcal{I}_K$, of size as equal as possible
\FOR{fold $k=1,\ldots,K$}
\STATE Label the covariates in batch $t$, fold $k$ by $X_{t1}^{(k)},\ldots,X_{tn_{t,k}}^{(k)}$
\IF{$t=1$}
\STATE Set $(\hat{e}_{t1}^{(k)},\ldots,\hat{e}_{tn_{t,k}}^{(k)})=\big(e_1(X_{t1}^{(k)}),\ldots,e_1(X_{tn_{t,k}}^{(k)})\big)$
\ELSE
\STATE Compute $(\hat{e}_{t1}^{(k)},\ldots,\hat{e}_{tn_{t,k}}^{(k)})$ using right-hand side of~\eqref{eq:e_t_hat_aipw} resp.~\eqref{eq:e_t_hat_epl} given batch $t$ budget constraints $m_{L,t} \leqslant m_{H,t}$
\ENDIF
\STATE Draw $(U_{t1}^{(k)},\ldots,U_{tn_{t,k}}^{(k)}) \stackrel{\mathrm{iid}}{\sim} \textnormal{Unif}(0,1)$
\STATE Assign treatments according to $Z_{ti}^{(k)} = \bm{1}(U_{ti}^{(k)} \leqslant \hat{e}_{ti}^{(k)})$, $i=1,\ldots,n_{t,k}$
\STATE Choose and store any $\hat{e}_t^{(k)}(\cdot) \in \mathcal{E}$ with $\hat{e}_t^{(k)}(X_{ti}^{(k)})=\hat{e}_{ti}^{(k)}, i=1,\ldots,n_{t,k}$
\ENDFOR
\STATE Observe outcomes $\{W_{ti} \mid 1 \leqslant i \leqslant N_t\}$ in batch $t$
\STATE Compute and store $N^{1/4}$-consistent estimates $\big(\hat{v}_{1:t}^{(1)}(\cdot,\cdot),\ldots,\hat{v}_{1:t}^{(K)}(\cdot,\cdot)\big)$ of $v_0(\cdot,\cdot)$ where for each fold $k=1,\ldots,K$,
$\hat{v}_{1:t}^{(k)}(\cdot,\cdot)$ depends only on the observations $\{W_{ui}:(u,i) \in \mathcal{I}_k, 1 \leqslant u \leqslant t\}$ in fold $k$ and batches $1,\ldots,t$
\ENDFOR
\STATE Compute final estimator $\hat\theta$ via~\eqref{eq:theta_hat}, using $(s_{a,\textnormal{AIPW}},s_{b,\textnormal{AIPW}})$, resp. $(s_{a,\textnormal{PL}},s_{b,\textnormal{PL}})$.
\end{algorithmic}
\end{algorithm}
For simplicity,
we assume treatment in the first batch is assigned according to a non-random propensity
$e_1(\cdot) \in \mathcal{F}_{\epsilon_1}$ for some $\epsilon_1>0$.
We let $\hat{e}_1^{(k)}(\cdot)=e_1(\cdot), k=1,\ldots,K$.
Then for later batches $t=2,\ldots,T$,
the target propensities are taken to be one of the following,
for an information function $\Psi(\cdot)$ satisfying Assumption~\ref{assump:Psi}:
\begin{equation}
\label{eq:e_t_star_aipw}
e_{t,\textnormal{AIPW}}^*(\cdot) \in \operatorname*{arg\,max}_{e_t(\cdot) \in \mathcal{F}_{*,t}} \Psi(V_{0:t,\textnormal{AIPW}}^{-1})
\quad\text{or}\quad
e_{t,\textnormal{EPL}}^*(\cdot) \in \operatorname*{arg\,max}_{e_t(\cdot) \in \mathcal{F}_{*,t}} \Psi(V_{0:t,\textnormal{EPL}}^{-1}).
\end{equation}
Above,
$V_{0:t,\textnormal{AIPW}}$ and $V_{0:t,\textnormal{EPL}}$ are the asymptotic variances of the oracle pooled estimators $\hat{\theta}_{\textnormal{AIPW}}$ and $\hat{\theta}_{\textnormal{EPL}}$ of~\eqref{eq:theta_star_aipw} and~\eqref{eq:theta_star_epl},
respectively,
\emph{when computed using observations in a non-adaptive batch experiment with only $t$ batches and propensities $e_1(\cdot),\ldots,e_t(\cdot)$}.
By~\eqref{eq:v_0_aipw} we can compute
\begin{equation}
\label{eq:V_0_t_aipw}
V_{0:t,\textnormal{AIPW}} = V_{0:t,\textnormal{AIPW}}(e_t;\eta_{0,\textnormal{AIPW}}) = \mathbb{E}_{P^X}\bigg[\frac{v_0(1,X)}{e_{0:t}(X)} + \frac{v_0(0,X)}{1-e_{0:t}(X)} + (\tau_0(X)-\theta_0)^2\bigg]
\end{equation}
where $\eta_{0,\textnormal{AIPW}}(x)$ includes the components $(v_0(0,x),v_0(1,x),\tau_0(x),\theta_0)$.
Similarly,
by~\eqref{eq:v_0_pl} we have
\begin{equation}
\label{eq:V_0_t_pl}
V_{0:t,\textnormal{EPL}}(e_t;\eta_{0,\textnormal{EPL}}) = \Bigg(\mathbb{E}_{P^X}\bigg[\frac{e_{0:t}(X)(1-e_{0:t}(X))}{v_0(0,X)e_{0:t}(X)+v_0(1,X)(1-e_{0:t}(X))}\psi(X)\psi(X)^{\top}\bigg]\Bigg)^{-1}
\end{equation}
where $\eta_{0,\textnormal{EPL}}(x)$ includes the components $(v_0(0,x),v_0(1,x))$.
In both of the preceding equations,
the dependence on the batch $t$ propensity score $e_t(\cdot)$ is through the mixture $e_{0:t}(\cdot)$ given by
\begin{equation}
\label{eq:e_0_x}
e_{0:t}(x) := \left(\sum_{u=1}^t \kappa_u\right)^{-1} \left(\sum_{u=1}^{t-1} \kappa_ue_u(x) + \kappa_te_t(x)\right), \quad x \in \mathcal{X}.
\end{equation}
Finally, the optimization set $\mathcal{F}_{*,t}$ in~\eqref{eq:e_t_star_aipw} is
\begin{equation}
\label{eq:F_t}
\mathcal{F}_{*,t}=\mathcal{F}_*(m_{L,t},m_{H,t};P^X)=\{e(\cdot) \in \mathcal{E} \mid m_{L,t} \leqslant \mathbb{E}_{P^X}[e(X)] \leqslant m_{H,t}\}
\end{equation}
which satisfies all the conditions in Assumption~\ref{assump:info_matrix} with covariate distribution $P^X$,
budget constraints $(m_{L,t},m_{H,t})$,
and base propensity class $\mathcal{E}$.
We do not target the final covariances $V_{0,\textnormal{AIPW}}=V_{0:T,\textnormal{AIPW}}$ and $V_{0,\textnormal{EPL}}=V_{0:T,\textnormal{EPL}}$ in our discussion here
since the sample sizes and budget constraints in future batches may not be known.
If they are known in advance,
then we can learn propensities for all future batches simultaneously at the time batch 2 covariates are observed,
and indeed target $V_{0:T,\textnormal{AIPW}}$ or $V_{0:T,\textnormal{EPL}}$ at that stage.
Now suppose we split our observations into $K$ folds as in a CSBAE,
and for notational simplicity we re-index the covariates in each batch $t=1,\ldots,T$, fold $k=1,\ldots,K$ as $X_{t1}^{(k)},\ldots,X_{tn_{t,k}}^{(k)}$.
Then following~\eqref{eq:e_hat},
we can estimate $e_{t,\textnormal{AIPW}}^*(\cdot)$ for each batch $t \geqslant 2$
within
each fold $k=1,\ldots,K$
by computing
\begin{equation}
\label{eq:e_t_hat_aipw}
\Big(\hat{e}_{t1,\textnormal{AIPW}}^{(k)},\ldots,\hat{e}_{tn_{t,k},\textnormal{AIPW}}^{(k)}\Big) \in \operatorname*{arg\,max}_{(e_1,\ldots,e_{n_{t,k}}) \in F_{n_{t,k}}} \Psi\biggl(\Big(\hat{V}_{0:t,\textnormal{AIPW}}^{(k)}\Big)^{-1}\bigg), \quad k=1,\ldots,K.
\end{equation}
Here $F_{n_{t,k}}=F_{n_{t,k}}(m_{L,t},m_{H,t})$ is defined as in~\eqref{eq:F_n},
and the estimate $\hat{V}_{0:t,\textnormal{AIPW}}^{(k)}$ of $V_{0:t,\textnormal{AIPW}}$ is given by
\begin{align}
\hat{V}_{0:t,\textnormal{AIPW}}^{(k)} & = \hat{V}_{0:t,\textnormal{AIPW}}^{(k)}\Bigl(e_1,\ldots,e_{n_{t,k}};\hat{e}_1^{(k)}(\cdot),\ldots,\hat{e}_{t-1}^{(k)}(\cdot),\hat{v}_{1:(t-1)}^{(k)}(\cdot,\cdot)\Bigr) \nonumber \\
& = -\frac{1}{n_{t,k}} \sum_{i=1}^{n_{t,k}} \frac{\hat{v}_{1:(t-1)}^{(k)}\big(1,X_{ti}^{(k)}\big)}{\hat{e}_{(0:t)i}^{(k)}} + \frac{\hat{v}_{1:(t-1)}^{(k)}\big(0,X_{ti}^{(k)}\big)}{1-\hat{e}_{(0:t)i}^{(t)}} \label{eq:V_0_t_hat_ate}
\end{align}
where each estimate $\hat{v}_{1:(t-1)}^{(k)}(z,\cdot)$ of the variance function $v_0(z,\cdot)$
is computed using only the observations from batches $u=1,\ldots,t-1$ within fold $k$.
Note that any plug-in estimate of $\tau_0(\cdot)$ and $\theta_0$ does not affect the optimization~\eqref{eq:e_t_hat_aipw}
and so can be omitted.
The dependence of $\hat{V}_{0:t,\textnormal{AIPW}}^{(k)}$ on the optimization variables $e_1,\ldots,e_{n_{t,k}}$ is through the mixture
quantities
\begin{equation}
\label{eq:e_0_i_hat}
\hat{e}_{(0:t)i}^{(k)} := \frac1{N_{1:t}} \Bigg(\,\sum_{u=1}^{t-1} N_u\hat{e}_u^{(k)}(X_{ti}^{(k)}) + N_te_i\Bigg), \quad i = 1,\ldots,n_{t,k}
\end{equation}
where for $u=1,\ldots,t-1$,
$\hat{e}_u^{(k)}(\cdot)$ is the (possibly adaptive) propensity used to assign treatment in batch $u$, fold $k$.
By comparing~\eqref{eq:e_0_i_hat} and~\eqref{eq:e_0_x},
we see that for each batch $u=1,\ldots,t$,
$N_u/N$ is being used as a plug-in estimate of $\kappa_u$
and the adaptive propensity score $\hat{e}_u^{(k)}(\cdot)$ is used as a plug-in estimate of its limit $e_u(\cdot)$.
Finally,
as in the conclusion of Lemma~\ref{lemma:concave_maximization},
the learned adaptive propensity $\hat{e}_{t,\textnormal{AIPW}}^{(k)}(\cdot)$ to be used for treatment assignment in batch $t$, fold $k$
is taken to be any choice in the predetermined base collection $\mathcal{E}$ with $\hat{e}_{t,\textnormal{AIPW}}^{(k)}\left(X_{ti}^{(k)}\right)=\hat{e}_{ti,\textnormal{AIPW}}^{(k)}$ for all $i=1,\ldots,n_{t,k}$.
Learning $e_{t,\textnormal{EPL}}^*(\cdot)$ is exactly analogous;
first we compute
\begin{equation}
\label{eq:e_t_hat_epl}
\Big(\hat{e}_{t1,\textnormal{EPL}}^{(k)},\ldots,\hat{e}_{tn_{t,k},\textnormal{EPL}}^{(k)}\Big) \in \operatorname*{arg\,max}_{(e_1,\ldots,e_{n_{t,k}}) \in F_{n_{t,k}}} \Psi\bigg(\Big(\hat{V}_{0:t,\textnormal{EPL}}^{(k)}\Big)^{-1}\bigg)
\end{equation}
where
\begin{align}
\label{eq:V_0_t_hat_pl}
\hat{V}_{0:t,\textnormal{EPL}}^{(k)} & =
\hat{V}_{0:t,\textnormal{EPL}}^{(k)}\Bigl(e_1,\ldots,e_{n_{t,k}};\hat{e}_1^{(k)}(\cdot),\ldots,\hat{e}_{t-1}^{(k)}(\cdot),\hat{v}_{1:(t-1)}^{(k)}(\cdot,\cdot)\Bigr) \\
& = \frac{1}{n_{t,k}} \sum_{i=1}^{n_{t,k}}\frac{\hat{e}_{(0:t)i}^{(k)}\big(1-\hat{e}_{(0:t)i}^{(k)}\big)}{\hat{v}_{1:(t-1)}^{(k)}\big(0,X_{ti}^{(k)}\big)\hat{e}_{(0:t)i}^{(k)} + \hat{v}_{1:(t-1)}^{(k)}\big(1,X_{ti}^{(k)}\big)\big(1-\hat{e}_{(0:t)i}^{(k)}\big)}\psi\big(X_{ti}^{(k)}\big)\psi\big(X_{ti}^{(k)}\big)^{\top}.
\end{align}
Then we assign treatment with any propensity $\hat{e}_{t,\textnormal{EPL}}^{(k)}(\cdot)$ in the base propensity class $\mathcal{E}$ satisfying
$\hat{e}_{t,\textnormal{EPL}}\big(X_{ti}^{(k)}\big)=\hat{e}_{ti,\textnormal{EPL}}^{(k)}$ for all $i=1,\ldots,n_{t,k}$.
The main additional regularity condition required to ensure the adaptive propensities $\hat{e}_{t,\textnormal{AIPW}}^{(k)}(\cdot)$ and $\hat{e}_{t,\textnormal{EPL}}^{(k)}$ above converge at the desired $O_p(N^{-1/4})$ RMS rate to $e_{t,\textnormal{AIPW}}^*$ and $e_{t,\textnormal{EPL}}^*$ of~\eqref{eq:e_t_star_aipw}
is the same rate of convergence in the estimates $\hat{\eta}^{(k)}_{\textnormal{AIPW}}(\cdot)$ and $\hat{\eta}^{(k)}_{\textnormal{EPL}}(\cdot)$ of the nuisance parameters $\eta_{0,\textnormal{AIPW}}(\cdot)$ and $\eta_{0,\textnormal{EPL}}(\cdot)$.
We also strengthen the sample size asymptotics~\eqref{eq:prop_asymp_limit} by requiring
\begin{equation}
\label{eq:prop_asymp_conv}
\frac{N_t}{N} = \kappa_t + O(N^{-1/4}), \quad t=1,\ldots,T.
\end{equation}
\begin{theorem}[Convergence of Algorithm~\ref{alg:csbae}]
\label{thm:concave_maximization}
Suppose $T \geqslant 2$,
fix a batch $t \in \{2,\ldots,T\}$
and suppose Assumption~\ref{assump:DGP} holds along with~\eqref{eq:prop_asymp_conv}.
Further assume treatment in batches $1,\ldots,t-1$ is assigned according to a CSBAE
where the batch 1 propensities are $\hat{e}_1^{(1)}(\cdot)=\ldots=\hat{e}_1^{(K)}(\cdot)=e_1(\cdot) \in \mathcal{F}_{\epsilon_1}$ for some $\epsilon_1 >0$,
and that for each fold $k=1,\ldots,K$:
\begin{itemize}
\item There exists an estimator $\hat{v}^{(k)}(\cdot)$ of the variance function $v_0(\cdot)$ depending only on the observations $\{W_{ui}^{(k)} \mid 1 \leqslant u \leqslant t-1\}$ in batches $1,\ldots,t-1$ assigned to fold $k$,
such that $\|\hat{v}^{(k)}(z,\cdot)-v_0(z,\cdot)\|_{2,P^X} = O_p(N^{-1/4})$ for $z=0,1$.
\item There are universal constants $0<c<C<\infty$ for which
\[
c \leqslant \inf_{(z,x) \in \{0,1\} \times \mathcal{X}} \min\big(\hat{v}^{(k)}(z,x),v_0(z,x)\big) \leqslant \sup_{(z,x) \in \{0,1\} \times \mathcal{X}} \max\big(\hat{v}^{(k)}(z,x),v_0(z,x)\big) \leqslant C.
\]
\item The information function $\Psi(\cdot)$ satisfies Assumption~\ref{assump:Psi}.
\end{itemize}
Let $\mathcal{E} \subseteq \mathcal{F}_0$ be any base propensity class satisfying the conditions of Assumption~\ref{assump:info_matrix},
and define $P_{N,t}^{(k),X}$ to be the empirical distribution on the covariates $X_{t1}^{(k)},\ldots,X_{tn_{t,k}}^{(k)}$ in batch $t$, fold $k$.
Then for any budget constraints $0 \leqslant m_{L,t} \leqslant m_{H,t} \leqslant 1$ and each fold $k=1,\ldots,K$,
the following holds:
\begin{enumerate}
\item (Design for $\hat{\theta}_{\textnormal{AIPW}}$)
There exists a target propensity $e_{t,\textnormal{AIPW}}^*(\cdot) \in \mathcal{E}$ satisfying~\eqref{eq:e_t_star_aipw}
that is unique $P^X$-almost surely.
Additionally, there exists a solution $(\hat{e}_{t1,\textnormal{AIPW}}^{(k)},\ldots,\hat{e}_{tn_{t,k},\textnormal{AIPW}}^{(k)})$ to~\eqref{eq:e_t_hat_aipw};
any such solution has the property that any propensity $\hat{e}_t^{(k)}(\cdot) \in \mathcal{E}$ with $\hat{e}_t^{(k)}(X_{ti}^{(k)})=\hat{e}_{ti,\textnormal{AIPW}}^{(k)}$ for $i=1,\ldots,n_{t,k}$ satisfies $$\|\hat{e}_t^{(k)}-e_{t,\textnormal{AIPW}}^*\|_{2,P^X} + \|\hat{e}_t^{(k)}-e_{t,\textnormal{AIPW}}^*\|_{2,P_{N,t}^X} = O_p(N^{-1/4}).$$
\end{enumerate}
If additionally,
the linear treatment effect assumption~\eqref{eq:pl_assumption} holds for some basis function $\psi(X)$ containing an intercept,
then:
\begin{enumerate}[resume]
\item (Design for $\hat{\theta}_{\textnormal{EPL}}$) There exists a target propensity $e_{t,\textnormal{EPL}}^*(\cdot) \in \mathcal{E}$ satisfying~\eqref{eq:e_t_star_aipw}
that is unique $P^X$-almost surely.
Furthermore,
there exists a solution $(\hat{e}_{t1,\textnormal{EPL}}^{(k)},\ldots,\hat{e}_{tn_{t,k},\textnormal{EPL}}^{(k)})$ to~\eqref{eq:e_t_hat_epl};
any such solution has the property that any propensity $\hat{e}_t^{(k)}(\cdot) \in \mathcal{E}$ with $\hat{e}_t^{(k)}(X_{ti}^{(k)})=\hat{e}_{ti,\textnormal{EPL}}^{(k)}$ for $i=1,\ldots,n_{t,k}$ satisfies $$\|\hat{e}_t^{(k)}-e_{t,\textnormal{EPL}}^*\|_{2,P^X} + \|\hat{e}_t^{(k)}-e_{t,\textnormal{EPL}}^*\|_{2,P_{N,t}^X} = O_p(N^{-1/4}).$$
\end{enumerate}
\end{theorem}
\begin{proof}
See Appendix \ref{proof:thm:concave_maximization}.
\end{proof}
\section{Numerical simulations}
\label{sec:simulations}
We implement Algorithm~\ref{alg:csbae} to construct some synthetic CSBAE's that illustrate the finite sample performance of our proposed methods.
For simplicity we consider $T=2$ batches throughout.
Our evaluation metric is
the average mean squared error (AMSE) of the estimators $\hat{\theta}_{\textnormal{AIPW}}$ and $\hat{\theta}_{\textnormal{EPL}}$
computed at the end of each CSBAE.
As a baseline,
we also compute feasible variants of the linearly aggregated estimators $\hat{\theta}_{\textnormal{AIPW}}^{(\mathrm{LA})}$ and $\hat{\theta}_{\textnormal{EPL}}^{(\mathrm{LA})}$ computed on a ``simple RCT":
that is,
a non-adaptive batch experiment with a constant propensity score in each batch.
We additionally consider the approach of~\citet{hahn2011adaptive} for design and estimation of $\theta_{0,\textnormal{ATE}}$.
This is equivalent to using Algorithm~\ref{alg:csbae} for design (without sample splitting, i.e.\ $K=1$)
and using the pooled $\hat{\theta}_{\textnormal{AIPW}}$ as the final estimator,
but with the covariates $X$ replaced everywhere by a coarse discretization $S=S(X)$.
As a hybrid we also consider using the discretized covariates $S$ for design but computing the final estimates $\hat{\theta}_{\textnormal{AIPW}}$
and $\hat{\theta}_{\textnormal{EPL}}$
using the full original covariate $X$.
To separately attribute efficiency gains to design and pooling,
we also consider the linearly aggregated estimators $\hat{\theta}_{\textnormal{AIPW}}^{(\mathrm{LA})}$ and $\hat{\theta}_{\textnormal{EPL}}^{(\mathrm{LA})}$ when a modification of Algorithm~\ref{alg:csbae} that targets $V_{2,\textnormal{AIPW}}$ and $V_{2,\textnormal{EPL}}$
(the asymptotic variances of the estimators $\hat{\theta}_{2,\textnormal{AIPW}}$ and $\hat{\theta}_{2,\textnormal{EPL}}$ given in Section~\ref{sec:pooled_estimation},
depending only on observations in batch 2)
is used for design.
We consider four data generating processes (DGPs),
distinguished by whether the covariate dimension $d$ is 1 or 10 and whether the conditional variance functions $v_0(\cdot,\cdot)$ are homoskedastic (with $v_0(0,x)=v_0(1,x)=1$ for all $x \in \mathcal{X}$) or heteroskedastic (with $v_0(0,x)=v_0(1,x)/2=\exp((\bm{1}_d^{\top}x)/(2\sqrt{d}))$.
The scaling by the covariate dimension $d$ in the heteroskedastic variance functions ensures the variance in $v_0(z,X)$ is independent of the covariate dimension $d$.
In all of the DGPs the covariates are i.i.d. spherical Gaussian,
i.e.\ $P^X=\mathcal{N}(0,I_d)$.
The outcome mean functions are taken to be $m_0(0,x)=m_0(1,x)=\bm{1}_d^{\top}x$ for $\bm{1}_d=(1,\ldots,1)' \in \mathbb{R}^d$.
For estimating $\theta_{0,\textnormal{PL}}$
we use the basis functions $\psi(x)=(1,x^{\top})^{\top} \in \mathbb{R}^p$ where $p=d+1$.
Note $\theta_{0,\textnormal{ATE}}=\theta_{0,\textnormal{PL}}=0$.
The potential outcomes $Y(0),Y(1)$ are generated as follows:
\begin{align*}
(\epsilon(0),\epsilon(1)) \mid X & \sim \mathcal{N}((0,0)^\top,\textnormal{diag}(v_0(0,X),v_0(1,X))) \\
Y(z) & = m_0(z,X) + \epsilon(z), \quad z=0,1.
\end{align*}
For each DGP we run Algorithm~\ref{alg:csbae} with $K=2$ folds,
batch sample sizes $N_1=N_2=1000$,
information function $\Psi=\Psi_a(\cdot)$ (corresponding to $A$-optimality),
and treatment fraction constraints $m_{L,t}=m_{H,t}=0.2$ for $t=1,2$.
In Appendix~\ref{app:simulations},
we present additional simulation results where $m_{L,1}=m_{H,1}$ remain at 0.2 but $m_{L,2}=m_{H,2}=0.4$,
so that the treatment budget for the second batch has increased.
The initial propensity score $e_1(\cdot)$ is taken to be constant (i.e. $e_1(x)=0.2$ for all $x \in \mathcal{X}$),
and the base propensity class $\mathcal{E}$ is the set of all $1$-Lipschitz functions taking values in $[0,1]$ when $d=1$ (cf. Example~\ref{example:lipschitz}).
When $d=10$,
we take $\mathcal{E}$ as in~\eqref{eq:E_conv_hull},
with $\{\theta_1,\ldots,\theta_m\}$ the collection of vectors $(a_1,\ldots,a_{11})' \in \mathbb{R}^{11}$ with each coordinate $a_i \in \{-2,-1,0,1,2\}$ and no more than two of the $a_i$'s nonzero.
The discretization used to implement the approach of~\citet{hahn2011adaptive} partitions $\mathbb{R}^d$ into four bins based on the quartiles of $\bm{1}_d^{\top}X$.
Note that this partition is along
the single dimension along which the variance functions $v_0(\cdot,\cdot)$ vary in the heteroskedastic DGPs,
so we would expect this to perform better than in practice,
where the structure of the variance functions is not known
(it could possibly be learned,
as in~\citet{tabord-meehan2022stratification}).
The choice of four bins is based on the experiments of~\citet{hahn2011adaptive},
which find minimal performance difference between two and six bins for the DGP's they consider.
We denote their ``binned" AIPW estimator by $\hat{\theta}_{\textnormal{AIPW}}^{(\text{bin})}$.
Recall,
as indicated above,
that this is equivalent to $\hat{\theta}_{\textnormal{AIPW}}$ when the only available covariate is the binned $S(X)$ and no cross-fitting is used.
For $\hat{\theta}_{\textnormal{AIPW}}^{(\text{bin})}$ the conditional means $\tilde{m}_0(z,s) = \mathbb{E}(Y(z) \mid S=s), z=0,1$ are estimated nonparametrically by sample outcome means among all units with $Z=z$ and $S=s$ in the appropriate batch(es) and fold(s).
Similarly, the conditional variances $\tilde{v}_0(z,s)=\textnormal{Var}(Y(z) \mid S=s)$ are estimated by sample outcome variances.
For all simulations,
the concave maximizations are performed using the CVXR software~\citep{fu2017cvxr} and the MOSEK solver~\citep{mosek}.
All estimates $\hat{m}(z,\cdot)$ of the mean function $m_0(z,\cdot), z=0,1$ are computed
by fitting a generalized additive model (GAM) to the outcomes $Y$ and covariates $X$ from the observations with treatment indicators equal to $z$ in the appropriate fold(s) and batch(es).
The GAMs use a thin-plate regression spline basis~\citep{wood2003thin-plate}
and the degrees of freedom are chosen using the generalized cross-validation procedure implemented in the \texttt{mgcv} package in R~\citep{wood2004stable}.
Variance function estimates $\hat{v}(z,\cdot)$
are computed by first computing $\hat{m}(z,\cdot)$ as above on the appropriate observations,
then fitting a GAM on these same observations to predict the squared residuals $(Y-\hat{m}(z,X))^2$ from $X$.
\subsection{Average treatment effect}
\begin{table}[!htb]
\centering
\caption{Simulated and asymptotic relative efficiencies (as defined in Section~\ref{sec:simulations} of the various design and estimation approaches for $\theta_{0,\textnormal{ATE}}$ that we study numerically under each of the four DGP's described in the text.
}
\label{table:ate_sim}
\begin{tabular}{ccccc}
\toprule
DGP & Estimator & Design & Sim. rel. eff. (90\% CI) & Asymp. rel. eff. \\
\midrule
\multirow{5}{0.15\linewidth}{\centering $d=1$, Homoskedastic} & $\hat{\theta}_{\textnormal{AIPW}}$ & Flexible & 0.989 (0.965, 1.013) & 1.000 \\
& $\hat{\theta}_{\textnormal{AIPW}}$ & Binned & 1.011 (0.977, 1.046) & 0.999 \\
& $\hat{\theta}_{\textnormal{AIPW}}$ & Simple RCT & 1.016 (1.003, 1.029) & 1.000 \\
& $\hat{\theta}_{\textnormal{AIPW}}^{(\mathrm{LA})}$ & Flexible & 0.984 (0.971, 0.998) & 0.999 \\
& $\hat{\theta}_{\textnormal{AIPW}}^{(\text{bin})}$ & Binned & 0.919 (0.877, 0.961) & 0.885 \\
\midrule
\multirow{5}{0.15\linewidth}{\centering $d=1$, Heteroskedastic} & $\hat{\theta}_{\textnormal{AIPW}}$ & Flexible & 1.016 (0.968, 1.064) & 1.048 \\
& $\hat{\theta}_{\textnormal{AIPW}}$ & Binned & 1.046 (0.997, 1.096) & 1.041 \\
& $\hat{\theta}_{\textnormal{AIPW}}$ & Simple RCT & 0.997 (0.979, 1.016) & 1.000 \\
& $\hat{\theta}_{\textnormal{AIPW}}^{(\mathrm{LA})}$ & Flexible & 1.008 (0.971, 1.045) & 1.024 \\
& $\hat{\theta}_{\textnormal{AIPW}}^{(\text{bin})}$ & Binned & 1.001 (0.946, 1.057) & 0.963 \\
\midrule
\multirow{5}{0.15\linewidth}{\centering $d=10$, Homoskedastic} & $\hat{\theta}_{\textnormal{AIPW}}$ & Flexible & 1.039 (0.990, 1.089) & 1.000 \\
& $\hat{\theta}_{\textnormal{AIPW}}$ & Binned & 1.011 (0.959, 1.065) & 1.000 \\
& $\hat{\theta}_{\textnormal{AIPW}}$ & Simple RCT & 1.081 (1.052, 1.110) & 1.000 \\
& $\hat{\theta}_{\textnormal{AIPW}}^{(\mathrm{LA})}$ & Flexible & 1.047 (1.018, 1.077) & 1.000 \\
& $\hat{\theta}_{\textnormal{AIPW}}^{(\text{bin})}$ & Binned & 0.477 (0.439, 0.519) & 0.417 \\
\midrule
\multirow{5}{0.15\linewidth}{\centering $d=10$, Heteroskedastic} & $\hat{\theta}_{\textnormal{AIPW}}$ & Flexible & 1.081 (1.023, 1.143) & 1.036 \\
& $\hat{\theta}_{\textnormal{AIPW}}$ & Binned & 1.050 (0.986, 1.117) & 1.000 \\
& $\hat{\theta}_{\textnormal{AIPW}}$ & Simple RCT & 1.102 (1.070, 1.135) & 1.000 \\
& $\hat{\theta}_{\textnormal{AIPW}}^{(\mathrm{LA})}$ & Flexible & 0.991 (0.952, 1.030) & 1.024 \\
& $\hat{\theta}_{\textnormal{AIPW}}^{(\text{bin})}$ & Binned & 0.651 (0.597, 0.708) & 0.595 \\
\bottomrule
\end{tabular}
\end{table}
Table~\ref{table:ate_sim} shows the performance of the various design and estimation procedures for $\theta_{0,\textnormal{ATE}}$ that we consider.
Each entry is a ``relative efficiency":
that is, a ratio of the MSE of the baseline approach
(which computes the linearly aggregated $\hat{\theta}_{\textnormal{AIPW}}^{(\mathrm{LA})}$ on a simple RCT)
to the MSE of the relevant approach.
The simulated relative efficiencies in the table estimate these MSE's by averaging the squared error of each estimator over 1,000 simulations.
The 90\% confidence intervals for the true finite sample relative efficiency are computed using 10,000 bootstrap replications of these 1,000 simulations.
Finally, the asymptotic relative efficiencies in Table~\ref{table:ate_sim} are computed by estimating the asymptotic variance of each estimator using the appropriate formula,
i.e. $V_{0,\textnormal{AIPW}}$ for $\hat{\theta}_{\textnormal{AIPW}}$
and $V_{\textnormal{AIPW}}^{(\mathrm{LA})}$ for $\hat{\theta}_{\textnormal{AIPW}}^{(\mathrm{LA})}$.
The computation is based on a non-adaptive batch experiment with second batch propensity $e_2(\cdot)$ equal to the average of the learned propensities from the relevant approach across the 1,000 simulations,
as a closed form solution for the limiting $e_2(\cdot)$ is not easily obtained in general.
All expectations over the covariate distribution are computed using Monte Carlo integration.
For both of the homoskedastic DGPs,
it is straightforward to show using Jensen's inequality that $e_2^*(x)=0.2$ for all $x$.
It can be further shown that the unpooled and pooled estimators are asymptotically equivalent.
Thus,
for the homoskedastic DGPs,
there is no asymptotic efficiency gain to be had.
With unequal budget constraints
there is some asymptotic benefit to pooling with homoskedastic variance functions (Appendix~\ref{app:simulations_unequal});
however the optimal design will still be the simple RCT.
Nonetheless,
we do observe some finite sample benefits to pooling in Table~\ref{table:ate_sim}.
For example,
the simulated relative efficiency of $\hat{\theta}_{\textnormal{AIPW}}$ on the simple RCT
is significantly larger than 1 for both $d=1$ and $d=10$.
We attribute this to improved nuisance estimates for the pooled estimator,
as discussed further in Appendix~\ref{app:oracle_sim}.
As one might expect,
this finite sample improvement from pooling is apparently offset by variance in both the flexible and binned design procedures.
Still,
in the homoskedastic DGPs,
our adaptive approaches do not show significant finite sample performance decline relative to the baseline,
which here is an oracle.
With unequal treatment constraints,
we further obtain asymptotic efficiency gains from pooling (Appendix~\ref{app:simulations_unequal}).
We notice that using the discretized covariate $S(X)$ in place of $X$ for both design and estimation,
as in~\citet{hahn2011adaptive},
leads to a substantial loss of efficiency.
Indeed, when $d=10$ the (asymptotic and simulated) variance of the estimator $\hat{\theta}_{\textnormal{AIPW}}^{(\text{bin})}$ is more than double that of our baseline under the homoskedastic DGP,
both asymptotically and in our finite sample simulations.
This efficiency loss occurs because the discretized $S(X)$ explains much less of the variation in the potential outcomes $Y(z)$ than the original $X$.
We expect greater precision losses from this discretization at the estimation stage when the variance functions $v_0(z,\cdot)$ vary substantially within the strata defined by $S(X)$.
For the heteroskedastic DGPs,
we see modest asymptotic efficiency gains from both pooling and design.
Design using the flexible base propensity class leads to about a 2.4\% asymptotic efficiency gain for both $d=1$ and $d=10$,
while pooling provides an additional 1--2\% gain.
These small asymptotic gains appear to be largely canceled out at our sample sizes by the finite sample variability in learning propensity scores,
limiting the net finite sample gains from design.
Of course, with greater heteroskedasticity and/or differences between $v_0(0,\cdot)$ and $v_0(1,\cdot)$,
we would expect greater efficiency gains from design;
in our simulations we have chosen to keep these differences within common ranges in social science studies as per~\citet{blackwell2022batch}.
\subsection{Partially linear model}
Unlike for estimating $\theta_{0,\textnormal{ATE}}$,
for estimating $\theta_{0,\textnormal{PL}}$,
we see clear efficiency gains over the baseline from \emph{design} when $d=1$ (Table~\ref{table:epl_sim}).
For instance,
the linearly aggregated estimator $\hat{\theta}_{\textnormal{EPL}}^{(\mathrm{LA})}$ exhibits a 5.6\% asymptotic efficiency gain as a result of the flexible design;
replacing this with the pooled estimator $\hat{\theta}_{\textnormal{EPL}}$
then yields a total asymptotic gain of 10.0\% over the baseline,
even in the homoskedastic DGP.
The analogous gains for the heteroskedastic DGP are slightly larger.
We once again observe a substantial finite sample benefit to pooling,
with the simulated relative efficiency of the approaches using the pooled $\hat{\theta}_{\textnormal{EPL}}$ tending to be larger than the asymptotic relative efficiency.
We attribute this to both improved use of nuisance estimates by the pooled estimator
(as in the ATE case)
as well as a more fundamental finite sample efficiency boost due to the fact that the asymptotic variance $V_{0,\textnormal{EPL}}$ is not exact for the oracle pooled $\hat{\theta}_{\textnormal{EPL}}^*$ in finite samples,
whereas $V_{0,\textnormal{AIPW}}$ is exact for $\hat{\theta}_{\textnormal{AIPW}}^*$;
see Appendix~\ref{app:oracle_sim}.
When $d=10$,
design introduces some more salient finite sample variance from the errors in the concave maximization procedure.
This is offset in both the homoskedastic and heteroskedastic DGPs by the asymptotic gains from the flexible design,
so that when $\hat{\theta}_{\textnormal{EPL}}$ is used,
the flexible design ultimately performs similarly to the simple RCT in finite samples.
Even if we use $\hat{\theta}_{\textnormal{EPL}}$ and hence the original covariate(s) $X$ in the final estimation step,
we see that the ``binned" design which uses only $S(X)$ for choosing the second batch propensity score
struggles to learn a substantially better propensity score than the simple RCT in all DGPs.
Indeed, in the heteroskedastic DGP with $d=1$,
we see a 2.1\% asymptotic efficiency \emph{loss} relative to the baseline from using the binned design.
By comparison, there is an 11.2\% asymptotic efficiency gain from using the flexible design.
The efficiency loss can occur with the binned design because the objective in~\eqref{eq:e_star} changes when working in terms of the discretized covariate $S(X)$ instead of $X$.
In other words, even though the simple RCT propensity $e_2(x)=0.2$ is within the class of propensities that can be chosen by the binned design,
it is worse according to the binned objective based on $S(X)$,
but not according to the objective based on the original $X$. \\
\begin{table}[!htb]
\centering
\caption{Same as Table~\ref{table:ate_sim},
but for estimating $\theta_{0,\textnormal{EPL}}$.
\\
}
\label{table:epl_sim}
\begin{tabular}{ccccc}
\toprule
DGP & Estimator & Design & Sim. rel. eff. (90\% CI) & Asymp. rel. eff. \\
\midrule
\multirow{4}{0.15\linewidth}{\centering $d=1$, Homoskedastic} & $\hat{\theta}_{\textnormal{EPL}}$ & Flexible & 1.150 (1.103, 1.198) & 1.100 \\
& $\hat{\theta}_{\textnormal{EPL}}$ & Binned & 1.057 (1.009, 1.107) & 1.026 \\
& $\hat{\theta}_{\textnormal{EPL}}$ & Simple RCT & 1.038 (1.018, 1.058) & 1.000 \\
& $\hat{\theta}_{\textnormal{EPL}}^{(\mathrm{LA})}$ & Flexible & 1.049 (1.016, 1.083) & 1.056 \\
\midrule
\multirow{4}{0.15\linewidth}{\centering $d=1$, Heteroskedastic} & $\hat{\theta}_{\textnormal{EPL}}$ & Flexible & 1.310 (1.214, 1.412) & 1.112 \\
& $\hat{\theta}_{\textnormal{EPL}}$ & Binned & 1.049 (0.972, 1.133) & 0.979 \\
& $\hat{\theta}_{\textnormal{EPL}}$ & Simple RCT & 1.084 (1.014, 1.159) & 1.000 \\
& $\hat{\theta}_{\textnormal{EPL}}^{(\mathrm{LA})}$ & Flexible & 0.983 (0.880, 1.076) & 1.073 \\
\midrule
\multirow{4}{0.15\linewidth}{\centering $d=10$, Homoskedastic} & $\hat{\theta}_{\textnormal{EPL}}$ & Flexible & 1.229 (1.203, 1.256) & 1.023 \\
& $\hat{\theta}_{\textnormal{EPL}}$ & Binned & 1.177 (1.147, 1.206) & 1.000 \\
& $\hat{\theta}_{\textnormal{EPL}}$ & Simple RCT & 1.228 (1.204, 1.251) & 1.000 \\
& $\hat{\theta}_{\textnormal{EPL}}^{(\mathrm{LA})}$ & Flexible & 1.023 (1.011, 1.035) & 1.019 \\
\midrule
\multirow{4}{0.15\linewidth}{\centering $d=10$, Heteroskedastic} & $\hat{\theta}_{\textnormal{EPL}}$ & Flexible & 0.997 (0.963, 1.032) & 1.071 \\
& $\hat{\theta}_{\textnormal{EPL}}$ & Binned & 0.932 (0.899, 0.965) & 1.002 \\
& $\hat{\theta}_{\textnormal{EPL}}$ & Simple RCT & 0.982 (0.952, 1.014) & 1.000 \\
& $\hat{\theta}_{\textnormal{EPL}}^{(\mathrm{LA})}$ & Flexible & 0.969 (0.942, 0.996) & 1.060 \\
\bottomrule
\end{tabular}
\end{table}
\section{Discussion}
\label{sec:discussion}
We view our primary technical contribution in this paper to be a careful extension of the double machine learning framework that enables both estimation and design in batched experiments based on pooled treatment effect estimators.
This allows the investigator to take advantage of the efficiency gains from pooling and design without needing to make strong parametric assumptions or to discretize their covariates.
As our numerical study in Section~\ref{sec:simulations} shows,
the latter can more than wipe out any efficiency gains from design.
Related to our work is the extensive literature on combining observational data with a (single batch) randomized experiment.
In that setting a primary concern is mitigating bias from unobserved confounders in the observational data~\citep{rosenman2021designing,gagnon2023precise}.
By contrast,
in our setting unconfoundedness holds by design in each batch of the experiment.
It would be useful to examine if the ideas from the present work can be extended to the observational setting where confounding bias is a concern. \\
\noindent
\textbf{Acknowledgments:} The authors thank Stefan Wager, Lihua Lei, and Kevin Guo for comments that improved the content of this paper.
{H.L. was partially supported by the Stanford Interdisciplinary Graduate Fellowship (SIGF).
This work was also supported by the NSF under grant DMS-2152780.