EconBase
← Back to paper

Optimal Experimental Design for Staggered Rollouts

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

123,267 characters · 18 sections · 52 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Optimal Experimental Design for Staggered Rollouts

\RUNAUTHOR \RUNTITLE{Optimal Experimental Design for Staggered Rollouts} \TITLE{Optimal Experimental Design for Staggered Rollouts{ \footnote{We thank seminar participants at Boston University, Columbia, Cornell, Cornell Tech, Emory, University of Florida, London Business School, National University of Singapore, Stanford, Toronto Rotman, University of Washington, UT Austin, Yale, Lyft Rideshare Labs, and participants at several conferences. We thank the editor, associate editor, and two anonymous referees for their insightful and helpful comments. Alphabetical author order other than the first author. Athey and Imbens were supported by the Office of Naval Research under grant N00014-19-1-2468. Bayati was supported by the National Science Foundation grant CMMI: 1554140. }}}

\ARTICLEAUTHORS{ \AUTHOR{Ruoxuan Xiong} \AFF{Department of Quantitative Theory and Methods, Emory University, \EMAIL{[email removed]}} \AUTHOR{Susan Athey} \AFF{Graduate School of Business, Stanford University, \EMAIL{[email removed]}} \AUTHOR{Mohsen Bayati} \AFF{Graduate School of Business, Stanford University, \EMAIL{[email removed]}} \AUTHOR{Guido Imbens} \AFF{Graduate School of Business, Stanford University, \EMAIL{[email removed]}} }

\ABSTRACT{ In this paper, we study the design and analysis of experiments conducted on a set of units over multiple time periods where the starting time of the treatment may vary by unit. The design problem involves selecting an initial treatment time for each unit in order to most precisely estimate both the instantaneous and cumulative effects of the treatment. We first consider non-adaptive experiments, where all treatment assignment decisions are made prior to the start of the experiment. For this case, we show that the optimization problem is generally NP-hard and we propose a near-optimal solution. Under this solution the fraction entering treatment each period is initially low, then high, and finally low again. Next, we study an adaptive experimental design problem, where both the decision to continue the experiment and treatment assignment decisions are updated after each period's data is collected. For the adaptive case we propose a new algorithm, the Precision-Guided Adaptive Experiment (PGAE) algorithm, that addresses the challenges at both the design stage and at the stage of estimating treatment effects, ensuring valid post-experiment inference accounting for the adaptive nature of the design. Using realistic settings, we demonstrate that our proposed solutions can reduce the opportunity cost of the experiments by over 50%, compared to static design benchmarks.

}

\KEYWORDS{Adaptive Experiments, Treatment Effect Estimation, Cumulative Effects, Panel Data, Dynamic Programming}

\setcounter{equation}{0} {\thesection.\arabic{equation}} \oldsection{Introduction} Large technology companies run tens of thousands of experiments (also known as A/B tests) per year to evaluate the impact of various decisions, new features, or products gupta2019top. In many cases, outcomes are observed for multiple periods, where the units can change treatment status over time. We refer to these experiments as \cmtfinal{panel experiments}; variants of these are sometimes referred to as experiments with staggered rollout or staggered adoption athey2018matrix, or stepped wedge designs hemming2015stepped.

example[Driver Experience] Consider a ride-hailing platform that plans to test the impact of a new app feature that improves driver experience. They wish to run an experiment to estimate the effect of providing the app feature to all drivers. To avoid biases from interference between drivers it is useful to randomize at the city level where the starting time for the intervention may vary by city.
example[Public Health Intervention] Consider a country that aims to measure the effect of a new public health intervention ({\it e.g.}, encouraging the use of masks or social distancing policies) on the spread of an infectious disease abaluck2021impact. To account for spillovers experimentation should be performed at an aggregate level ({\it e.g.}, cities). To facilitate the estimation of cumulative effects it is useful to vary the starting date.

The primary objective of this paper is to propose experimental designs to optimize the precision of post-experiment estimates of instantaneous and \cmtfinal{lagged} effects of a treatment. To operationalize this, the experimenter commits in advance to an estimator which has a precision associated with each experimental design. The challenge is to find a design that optimizes this precision given the estimator.

We assume all units start in the control (no treatment) state at the initial time period. The design problem is to select, for each unit, the time period to begin treatment. We assume that units cannot switch back to control during the experiment, and thus remain treated once exposed to the treatment. This is a common setting.\footnote{The setting where units can arbitrarily switch between treatment and control is discussed in Sections (ref) and (ref), and turns out to be a simpler setting.} The treatment allocation time can vary across units, leading to a staggered treatment adoption or stepped wedge design.

Summary of Contributions

{\bf Non-adaptive experiments. } We first study the design of non-adaptive experiments, where both the number of units and time periods and treatment decisions are determined prior to the start of the experiment. We focus on the properties of the generalized least squares (GLS) estimator. We consider a general version of feasible GLS to estimate instantaneous and lagged treatment effects\footnote{Lagged treatment effect measures the effect of this period's treatment on a future period's outcome.} from observed outcomes that can be non-stationary. \cmtrx{The estimator determines the precision of estimated instantaneous and lagged effects. A linear combination of these precisions comprises the objective function for the experimental design problem. } Finding the optimal solution is generally NP-hard. We provide the analytical optimality conditions for the design, and propose an algorithm to choose a treatment design based on the optimality conditions. The precision of our selected design approximates the optimal objective (best achievable precision) within a multiplicative factor of $1+O(1/N^2)$ where $N$ is the number of units.

Our solution to the design problem of non-adaptive experiments has two prominent features. First, the fraction of treated units per time period takes an $S$-shaped curve; the treatment rolls out to units more slowly (over time) in the beginning and at the end, but rapidly in the middle, of the experiment. Second, the optimal design imposes this rollout pattern for each stratum, where strata are defined by groups with the same observed and estimated latent covariate values.

{\bf Adaptive experiments.} Next we study the design of adaptive experiments, where the number of units is fixed, but the experiment can be terminated early and treatment assignment decisions can be adaptively made after each period's data is collected. These experiments are useful when the pre-set duration is more than needed to attain a target precision of treatment effect estimates, and when treatment decisions cannot be optimally made ex ante. In our main contribution, we propose a new algorithm, the Precision-Guided Adaptive Experiment or PGAE algorithm, that adaptively terminates the experiment based on the estimated precision of the estimated treatment effect from partially observed results of the experiment. It employs dynamic programming \cmtfinal{to adaptively optimize} \cmtfinal{the speed of treatment rollout in} subsequent time periods. The resulting adaptive experiment achieves a target precision, using a shorter duration, or equivalently incurring a lower cost, than the non-adaptive experiment.

The adaptive nature of the experiment creates challenges arising from the fact that the outcomes and assignments that occur early in the experiment affect the treatment assignments to units later in the experiment. Thus, the treatment assignments are not independent of observed outcomes. We propose an estimation method based on sample splitting that ensures that the estimates obtained by PGAE are consistent and asymptotically converge to a normal distribution. An appealing property of the proposed estimation scheme we propose is that the final treatment effect estimation uses all of the data, incurring no efficiency loss compared to an oracle who would have access to the same design at the beginning of the experiment.

Finally, we illustrate the superior performance of our solutions, as compared to benchmarks, for non-adaptive and adaptive experiments through synthetic experiments based on four real data sets about flu occurrence rates, home medical visits, grocery expenditure, and Lending Club loans.

Related Literature

Our \cmtfinal{staggered rollout designs of panel experiments} are related to a number of other experimental designs. They are most closely related to the stepped wedge designs in clinical trials brown2006stepped. Stepped wedge designs sequentially roll out an intervention to clusters over several periods. Prior work on optimal stepped wedge design hussey2007design,hemming2015stepped,li2018optimal studies optimal treatment assignments of clusters under a linear mixed outcome model that has no observed or latent covariates and assumes constant treatment effect over time with instantaneous effects only. In contrast, we study the optimal design allowing treatment effects to vary over time ({\it i.e.}, with instantaneous and lagged effects).

Related to our design, a few other designs proposed recently are suitable for studying time-varying treatment effects. One design is the synthetic control design \cmtfinal{for} panel experiments \cmtfinal{where the design} selects \cmtfinal{units to be treated}, allocates treatment to \cmtfinal{all of} them \cmtfinal{in a single period}, and forms a synthetic treated and control unit for treatment effect estimation doudchenkodesigning2021,doudchenko2021synthetic,abadie2021synthetic. Another design is the switchback design bojinov2020design,xiong2023bias for a single experimental unit that can arbitrarily switch between treatment and control. In contrast to these two designs, our designs leverage variation in treatment times across units to increase power. The randomized designs proposed in basse2019minimax also allow for cross-unit variation in treatment times, but they are studied under a different framework that minimizes the worst-case risk in randomization-based inference.

When interference between units is a concern, our design can be modified by following a conventional approach to avoid biases by aggregating units to a level that \cmtfinal{interference} is no longer a problem. But we note a growing literature that directly tackles the biases using novel experimental design ideas, such as multiple randomization bajari2021multiple,johari2020experimental and designs \cmtfinal{that perturb treatments near equilibrium outcomes} wager2021experimenting. In contrast, our design, by abstracting away from interference, can be used to study the rich dynamics of cumulative effects over time.

Different from all the aforementioned designs, we additionally study the design and analysis of adaptive experiments. Our proposed PGAE algorithm consists of three components: adaptive treatment decisions, an experiment termination rule, and post-experiment inference. The adaptive treatment decision \cmtfinal{component} relates to the literature on adaptive designs in sequential experiments ({\it e.g.,} efron1971forcing,bhat2019near,glynn2020adaptive), online learning and multi-armed bandits ({\it e.g.,} bubeck2012regret,lattimore2018bandit). \cmtfinal{Distinct} from this literature, we consider a \cmtfinal{panel setting}, and our design choices allow for experiment termination and ensure that inference is manageable.

The experiment termination rule \cmtfinal{component} is based on the precision of the treatment effect estimate. \cmtrx{The precision-based termination rule has been used to obtain fixed-volume confidence sets from a sequence of independent random variables {glynn1992asymptotic,glynn1992asymptoticstopping,singham2012finite}. We extend the use of this rule to the panel setting.} Our proposed PGAE algorithm decides whether to terminate the experiment at every period, which is related to the sequential testing problem, considered by siegmund1985sequential,wald2004sequential,bertsekas2012dynamic,johari2017peeking,ju2019sequential among others. The key challenge in sequential testing is to draw valid inference post-experiment. The PGAE algorithm addresses this challenge through a sample-splitting technique.

We split units into disjoint sets, with each serving a different purpose. The idea of sample splitting has been used in the econometrics literature \cmtrx{(angrist1995split,angrist1999jackknife,athey2016recursive,chernozhukov2018double among others)} for valid inference when machine learning methods are used and overfitting is a concern. \cmtmb{In contrast, we use the full sample for the final treatment effect estimation that incurs no loss in estimation efficiency}. However, splitting samples into disjoint sets of units is crucial in decoupling experiment termination from inference in PGAE. In multi-armed bandits and other sequential decision problems, the splitting of \cmtfinal{data in adaptive experiments} has been considered by, among others, \cmtmb{auer2003using,goldenshluger2013linear,bastani2020online,hamidi2019personalizing} for consistent estimation of the arm parameters. In contrast to this literature, we repeatedly experiment on the same set of units, and split samples in the unit dimension. \cmtmb{Finally, we note that lai1982least showed normal approximation for the estimation error of adaptive least squares under a certain stability condition on the inverse covariance matrix. Recently, stronger results were obtained by leveraging debiasing techniques for OLS deshpande2017accurate and LASSO deshpande2021online. In our panel setting, the stability condition of lai1982least does not hold; however, the multi-unit nature of the problem allows us to avoid debiasing, even without the stability condition. Moreover, in our setting the experiment stopping rule is adaptively selected. }

\setcounter{equation}{0} {\thesection.\arabic{equation}} \oldsection{\cmtfinal{Panel Experiments with Staggered Rollouts} and Assumptions} In this section, we first introduce \cmtfinal{panel experiments with staggered rollouts} and then describe the two types of experiments that are considered in this paper. We also specify assumptions on the outcome model and treatment designs. Throughout, $[M]$ refers to the set $\{1, 2, \ldots, M\}$, for any positive integer $M$, {and additional notations used in the paper are summarized in Table (ref)}.

Recall Examples (ref)-(ref), and assume we are planning to run an experiment to estimate the effect of a treatment of interest (app feature or public health intervention). Let $z_{it}$ in $\{-1,+1\}$ be the treatment variable for unit $i$ at time $t$, for $i, t \in \mathbb{Z}$, where $z_{it} = +1$ means unit $i$ is treated at time $t$ and $z_{it} = -1$ means otherwise. Assume the treatment is not applied to any unit before the experiment starts, that is, $z_{it} \equiv -1$ for all $i$ and $t \leq 0$. The experimental designer decides the treatment assignments in an experiment with $N$ units over $T$ time periods, that is, chooses $z_{it}$ for $i \in [N]$ and $t \in [T]$. All the decision variables can be written in a compact form $Z = [z_{it}]_{(i,t)\in [N] \times [T]}$, which is referred to as the treatment design of the experiment. Different treatment designs lead to different panels of observed outcomes, denoted as $Y = [Y_{it}]_{(i,t)\in [N] \times [T]}$, that affect the precision of treatment effect estimates. The estimation precision, to be defined formally in Section (ref), directly impacts the cost of running the experiment since it specifies the required number of units and time periods. Therefore, optimizing the treatment decisions is an integral part of designing \cmtfinal{panel experiments}.

In this paper, we consider two types of \cmtfinal{panel experiments}: non-adaptive experiments and adaptive experiments. For nonadaptive experiments, $N$ and $T$ are fixed, whereas for adaptive experiments, $N$ is fixed, but $T$ is unknown. Non-adaptive experiments enjoy the benefit of simplicity and involve the straightforward construction of statistical tests for the treatment effects. However, non-adaptive experiments can be inefficient in time and cost, as the experiment may run longer than needed to attain a certain precision of treatment effect estimates. In comparison, adaptive experiments can early stop the experiment if needed, at the expense of making statistical inference for treatment effects more challenging. We study the design of non-adaptive experiments in Section (ref), where treatment decisions are made before the experiment starts. Building on our insights from the non-adaptive experiments, we propose an algorithm in Section (ref) to run adaptive experiments that make treatment decisions adaptively during the experiment and provide valid post-experiment inference of treatment effects.

\cmtfinal{We focus on a class of panel experiments whose treatment assignments satisfy the irreversible treatment adoption condition. These experiments with treatment adoption times possibly varying across units are referred to as panel experiments with staggered rollouts; the design of these experiments is referred to as the staggered rollout design.}

assumption[Irreversible Treatment Adoption] For all $(i,t)\in [N]\times[T]$, $z_{it} \leq z_{i,t+1}$.

We mostly focus on the irreversible scenario for two reasons. First, often there are practical constraints restricting units from switching between control and treatment, such as the policies or programs implemented at the group level ({\it e.g.}, city or state level), or the features that are electronically rolled to the user interface ({\it e.g.}, ride-hailing app feature). Second, when the treated units switch back to the control, they may not return to the original control status. For example, drivers may develop different driving habits via the new feature. In fact, the irreversible pattern is common in observational studies ({\it e.g.}, card1994minimum,abadie2010synthetic). In the experiment design literature, the irreversible pattern appears in the stepped wedge designs of cluster randomized trials in public health brown2006stepped,hussey2007design,woertman2013stepped,hemming2015stepped, and in the synthetic control designs doudchenkodesigning2021,doudchenko2021synthetic,abadie2021synthetic.

remark[Extension to Reversible Treatments] In Section (ref), we study the design of non-adaptive panel experiments when the treatment can be stopped. We show that optimizing the treatment designs for this case follows from our results in Section (ref).
figure[figure omitted — 1,036 chars of source]

Next we introduce the potential outcome model.

assumption[Potential Outcome Model] The potential outcomes for unit $i$ at time $t$ can be written as \cmtrx{\[Y_{it}(z_{i,t-\ell}, \cdots, z_{i,t-1}, z_{it}) \] for a nonnegative integer $\ell$, where $\ell$ is known. }

Assumption (ref) requires that a unit's treatment $z_{it}$ only affects this unit's own outcomes (no cross-unit interference). It further requires that a unit's outcomes are not affected by this unit's future treatment assignments (no anticipation). For example, drivers do not increase their working hours in anticipation of the new app feature. No anticipation is also commonly imposed in the literature basse2019minimax,bojinov2020design.

We allow the potential outcomes to depend on the treatment assignments up to $\ell$ periods in the past. \cmtrx{We assume $\ell$ is known for the design and estimation purposes. Note that there exists a bias-variance tradeoff in the specification of $\ell$ in practice. When the specified $\ell$ is less than the true value, the estimator can be biased. When the specified $\ell$ is more than the true value, the estimator is unbiased, but can be inefficient. See Figure (ref) in Section (ref) for an example. }

We define the average instantaneous effect $\tau_0$ and lagged effects $\tau_1,\dots,\tau_\ell$ as \cmtrx{ \[ \tau_j \coloneqq \frac{1}{NT} \sum_{i,t} \frac{1}{2} \Big[ Y_{it}(-1, \cdots, -1, \underbrace{+1}_{z_{i,t-j}}, +1, \cdots, +1) - Y_{it}(-1, \cdots, -1, \underbrace{-1}_{z_{i,t-j}}, +1, \cdots, +1) \Big]\,, \]} for all $j\in \{0\} \cup [\ell]$. Here $1/2$ is used for notation simplicity in the specification (ref) below and remaining sections to account for the scaling of $z_{it}$ that takes value between $\{+1,-1\}$ as opposed to $\{+1,0\}$. Note that $\tau_j$ is the average, over individuals and time periods, of individual-time specific treatment effects that can be heterogeneous in both units and time periods.

Here $\tau_0, \tau_1, \cdots, \tau_\ell$ can take arbitrary values, and therefore allow for general dynamics of cumulative effects over time. To demonstrate this generality, we describe a few scenarios below.

itemize• The average effect accumulates over time, for example when $\tau_0,\ldots,\tau_\ell$ all have the same sign, as shown in Figure (ref). • The average effect attenuates over time, for example when $\tau_0$ and $\tau_1$ have a different sign than all of $\tau_2,\ldots,\tau_\ell$, as shown in Figure (ref). • The average effect is constant over time, for example when $\ell = 0$. • The average cumulative effect is zero, for example when $\sum_{j=0}^\ell \tau_j=0$. After $\ell$ periods, a unit's outcome is the same as the outcome of never being treated.

In this paper, we are interested in estimating $\bm{\tau} = (\tau_0, \tau_1, \cdots, \tau_\ell)$. Below we study how to choose $Z = [z_{it}]_{(i,t)\in [N]\times[T]}$ to accurately estimate these parameters.

\setcounter{equation}{0} {\thesection.\arabic{equation}} \oldsection{Non-Adaptive Experiments} In this section, we study the \cmtfinal{staggered rollout} design of non-adaptive experiments. In Section (ref), we present an approach for estimating $\tau_0, \tau_1, \cdots, \tau_\ell$. Next, in Section (ref), we define the precision of the proposed estimators. We formulate an optimization problem, where for a given estimator, the solution yields the set of treatment times for each unit that maximizes the precision of the estimator. Finally, we present our solution to the optimization problem in Section (ref).

Estimation of Instantaneous and Lagged Effects

The decision maker needs to consider two objectives when choosing an estimator for $\tau_0, \tau_1, \cdots, \tau_\ell$. First is the statistical properties of this estimator, such as the bias, variance, and mean-squared error (MSE). Second is the feasibility of optimizing units' treatment times based on the properties of this estimator. For example, one can study the optimal treatment assignments based on MSE of a simple estimator \cmtrx{such as the difference-in-means estimator}, but this simple estimator can have a large \cmtrx{variance and} MSE. Alternatively, one can use a sophisticated estimator with a small MSE \cmtrx{such as regularized matrix factorization estimator in Section (ref)}, but the functional form of MSE may be very complex, which can lead to an intractable optimization problem in the design phase.

In this paper, we seek to balance these two objectives. We study the optimization of treatment assignments, for a particular estimator, namely the Generalized Least Squares (GLS) estimator, given a particular model for the data generating process. The GLS estimator for $\tau_0, \tau_1, \cdots, \tau_\ell$ is based on the following specification

equation[equation omitted — 239 chars of source]

where $\alpha_i$ and $\beta_t$ are unobserved unit and time fixed effects, $\*X_i \in \+R^{d_x}$ are observed covariates of dimension $d_x$, $\bm{\theta}_t \in\+R^{d_x}$ are their (non-random) unobserved time-varying coefficients $\bm{\theta}_t$, $\*u_i \in \+R^{d_u}$ are latent (non-random) covariates of dimension $d_u$, $\*v_t \in \+R^{d_u}$ are (random) latent factors, and $\{\varepsilon_{it}\}_{it\in[N]\times[T]}$ are unobserved residuals. We allow $d_x = 0$, that is, (ref) does not have observed covariates. We also allow $d_u = 0$, that is, (ref) does not have latent covariates.

The specification with two-way ({\it i.e.}, unit and time) fixed effects and possibly with observed covariates has been widely used in observational studies to estimate treatment effects angrist2008mostly.\footnote{There are two reasons. First, $\alpha_i$ and $\beta_t$ can be arbitrarily correlated with the observed explanatory variables $(\*X_i, {z}_{it}, \cdots, {z}_{i,t-\ell})$, so they can capture the unobserved additive unit-specific and time-specific confounders that jointly affect ${z}_{it}$ and $Y_{it}$ angrist2008mostly. Second, $\alpha_i$ captures the mean effect of unobservables on unit $i$'s outcome over time, and $\beta_t$ captures the mean effect of unobservables on units' outcomes at time $t$. } Observed covariates $\*X_i$ in (ref) are time-invariant and are not affected by the treatment. Examples of such covariates are an individual's gender and race or attributes of a city.\footnote{The time-invariance assumption is commonly used in difference-in-differences estimators card1994minimum and synthetic control estimators abadie2010synthetic in observational studies.} Using $\alpha_i$, $\beta_t$ and $\*X_i$ controls for the heterogeneity in units and time periods, that can effectively reduce variance in treatment effect estimation. Note that (ref) has the interactive latent factor structure $\*u_i^\top \*v_t$ that can capture the multiplicative unobserved effects, and is therefore more general than the specification with additive effects $\alpha_i$ and $\beta_t$ only. \cmtrx{We assume $\*v_t$ is random, but $\*u_i$ is nonrandom for the design purpose. Such an assumption is also imposed in a different literature that estimates the latent factors on large panels bai2002determining,bai2003inferential. }

GLS estimates $\tau_0, \cdots, \tau_\ell$, $\alpha_i$, $\beta_t$ and $\bm{\theta}_t$ by minimizing the weighted sum of squared residuals $e_{it}$ in (ref) (see Section (ref) for more details). GLS has two nice properties \cmtfinal{when the optimal weights are chosen, following} the Gauss-Markov Theorem (see Lemma (ref) in Section (ref)) under the strict exogeneity assumption in Assumption (ref) below. First, GLS is the best linear unbiased estimator (BLUE) of $\tau_0, \cdots, \tau_\ell$, meaning that this estimator has the smallest variance among all the linear unbiased estimators. Second, there is an explicit formula for the variance and precision of GLS in terms of $z_{it}$, which makes the treatment design problem tractable.

assumption[Error Structure] \cmtrx{$\*v_t \in \+R^{d_u}$ is i.i.d. in $t$} with $\+E[\*v_t \mid \bm{\xi}_{it}] = \bm{0}$ and $\+E[\*v_t \*v_s^\top \mid \bm{\xi}_{it}, \bm{\xi}_{js} ] = \bm\Sigma_{v} \cdot \bm{1}_{t = s}$ for all $i, j, t$ and $s$, where $\bm{\xi}_{it} = (\*X_i, z_{it}, \cdots, z_{i-\ell, t})$. Moreover, $\varepsilon_{it} $ is i.i.d. in $i$ and $t$ with $\+E[\varepsilon_{it} \mid \bm{\xi}_{it}] =0$ and $\+E[\varepsilon_{it} \varepsilon_{js} \mid \bm{\xi}_{it}, \bm{\xi}_{js} ] = \sigma_\varepsilon^2 \cdot \bm{1}_{i = j, t = s}$ for all $i, j, t$ and $s$. In addition, $\varepsilon_{it}$ is independent of $\*v_s$ for all $i, t$ and $s$.
remark[Machine learning heuristics for $\bm{\tau}$] Instead of focusing on the least squares estimator, one can look at other estimators, for example, machine learning estimators. We provide one such estimator in Section (ref) for $\bm{\tau}$ that does not rely on Assumption (ref) about $\*v_t$. This type of machine learning-based estimators are typically biased in the estimation of $\bm{\tau}$, but could have a smaller variance than the unbiased estimators. The optimization of $Z$ based on these estimators is generally not tractable, and therefore we do not emphasize them in this paper.
remark[Implication of Assumption (ref)] $\+E[\*v_t] = 0$ holds without loss of generality, as we can always project their mean to $\alpha_i$.\footnote{If $\+E[\*v_t] \neq 0$, let $\alpha^\dagger_i = \alpha_i + \*u_i^\top \+E[\*v_t] $ and $\*v^\dagger_t = \*v_t - \+E[\*v_t]$, and then $\*v^\dagger_t $ has mean zero. } Under Assumption (ref), $e_{it}$ satisfies \[\+E[e_{it}] = 0, \qquad \+E[e_{it}^2] = \*u_i^\top \*\Sigma_{v} \*u_i + \sigma_{\varepsilon}^2, \qquad \+E[e_{it} e_{jt}] = \*u_i^\top \*\Sigma_{v} \*u_j, \qquad \+E[e_{it}e_{js}] = 0, \quad \text{ for } \,\, t \neq s, \] implying that $e_{it}$ is correlated in the unit dimension, but it is uncorrelated in the time dimension.

\cmtfinal{

remark[Feasiblity of GLS] The weight matrix $\*W$ in GLS is optimal when it is proportional to the inverse covariance matrix of $\bm{e}_t = [e_{it}]_{i \in [N]}$, {\it i.e.}, $\big(\*U \*\Sigma_v \*U^\top + \sigma_\varepsilon^2 \*I_N \big)^{-1} $ with $\*U = [\*u_i]_{i \in [N]}$. As the optimal weight matrix is unknown, in practice, we can use feasible GLS, where we first use OLS, estimate the optimal weight matrix, and then use GLS with the estimated weight matrix.

}

remark[Treatment Effect Estimation with Heterogeneity] Specifications like (ref) have been commonly used to estimate average treatment effects in empirical studies in many domains, such as economics card1994minimum, operations ({\it e.g.}, cui2019learning,cachon2019does), and healthcare abaluck2021impact. Such specifications do not restrict treatment effects to be homogeneous. In Remark (ref) below, we discuss an alternative estimation approach when heterogeneous treatment effects are the objects of interest.

Optimization Problems for Treatment Decisions

In this section, we introduce the optimization problem to solve the optimal $Z$ based on the statistical properties of GLS. From a high level perspective, we are interested in precisely estimating $\bm{\tau}$. However, there are $\ell+1$ parameters in $\bm{\tau}$ and it is generally infeasible to find an $Z$ that simultaneously maximizes the precision of each of $\hat{\tau}_0, \cdots, \hat{\tau}_{\ell}$ for $\ell \geq 1$. Instead, one needs to consider an objective function that summarizes the precision of each of $\hat{\tau}_0, \cdots, \hat{\tau}_{\ell}$ into a scalar. There are two categories of objective functions that one may be interested in:

itemize• Balancing the precision/variance of each of $\hat{\tau}_0, \cdots, \hat{\tau}_{\ell}$ • Maximizing the precision of some linear combination of $\hat{\tau}_0, \cdots, \hat{\tau}_{\ell}$

Examples of the objective functions related to the first category include atkinson2007optimum

itemize• A-optimal design: minimizes the trace of $\mathrm{Var}(\hat{\bm{\tau}})$ • D-optimal design: minimizes the determinant of $\mathrm{Var}(\hat{\bm{\tau}})$ • T-optimal design: maximizes the trace of the inverse of $\mathrm{Var}(\hat{\bm{\tau}})$

An example of the objective functions related to the second category is

itemize• Cumulative effect ($\sum_{j=0}^\ell \tau_j$): minimizes $\mathrm{Var}\big(\sum_{j = 0}^\ell \hat{\tau}_j\big)$

When $\ell = 0$, the four objective functions mentioned above are equivalent to one another. For general $\ell$, they are not equivalent, and the optimal treatment assignments for these four objectives can be different (but the difference can be small).

In this section, we focus on finding the optimal assignment $A = [A_i]_{i \in [N]}$, where $A_i \in [T] \cup \{\infty\}$ denotes the first time that unit $i$ adopts the treatment ($A_i=\infty$ means unit $i$ was never treated\footnote{We use “$\infty$” to preserve the ordering that a larger value for $A_i$ implies the treatment is assigned at a later time.}). Because the treatment is irreversible, there is a one-to-one mapping between $(z_{i1}, z_{i2}, \cdots, z_{iT})$ and $A_i$. Given the reduction to the $N$ adoption times we focus on analytically solving the T-optimal (or trace-optimal) design:

equation[equation omitted — 111 chars of source]

\cmtrx{where $\mathrm{Prec}(\hat{\bm{\tau}}; A)$ is the precision matrix of $\bm{\tau}$ and is defined as the matrix inverse of the variance-covariance matrix of $\bm{\tau}$, $\mathrm{Var}(\hat{\bm{\tau}}; A)$, when the assignment $A$ is used.}

T-optimal design was first introduced by atkinson1975design to discriminate between two competing regression models ({\it e.g.}, to determine whether $\ell = \ell_0$ or $\ell = \ell_1$ is true, for distinct number of lags $\ell_0$ and $\ell_1$ in our setting). Since then, the T-optimal design has been studied by atkinson1975optimal,ucinski2005t,wiens2009robust,dette2012t,dette2013robust,dette2015bayesian and others.

Solving the T-optimal design is generally challenging, even in some special cases (see Example (ref)). We provide the explicit optimality conditions for the integer program (ref) in Section (ref). Based on the optimality conditions, we provide an algorithm on choosing a design in Algorithm (ref) in Section (ref). Admittedly, the other three objectives mentioned above would be natural and of practical interest, especially the one that minimizes the variance of the estimated cumulative effect, $\mathrm{Var}\big(\sum_{j = 0}^\ell \hat{\tau}_j\big)$. However, analytically solving the other three objectives is generally infeasible, as explained in Remark (ref) below. Instead, one could numerically solve the other three objectives in practice. We visualize the numerical solutions for D-optimal design in Figure (ref) in Section (ref), which has a similar structure as our solutions to (ref). We empirically show in Section (ref) and Section (ref) that our solutions to (ref) outperform the benchmark treatment designs measured by the objective of A-optimal design and $\mathrm{Var}\big(\sum_{j = 0}^\ell \hat{\tau}_j\big)$.

exampleIf $T=1$, $\ell = 0$, and all the covariates are observed ({\it i.e.}, $d_u = 0$), (ref) coincides with the offline optimization problem in bhat2019near, and is equivalent to the MAX-CUT problem and is NP-hard hayes2002computing,mertens2006easiest. In this paper, we focus on the case of low-dimensional covariates.\footnote{This makes sense as we allow for latent covariates, which can summarize the information and reduce the dimensionality of (high-dimensional) observed covariates.} Our solution from Algorithm (ref) is provably close to the optimal integer solution to (ref) for a large $N$.
remark[Challenges in alternative objectives] Note that each entry in $\mathrm{Var}(\hat{\bm{\tau}}) $ is a ratio of two polynomial functions of $z_{it}$, where the degrees of numerator and denominator are $2\ell$ and $2(\ell+1)$, respectively, for $\ell \geq 1$. The objective of A-optimal design and $\mathrm{Var}\big(\sum_{j = 0}^\ell \hat{\tau}_j\big)$ are both sums of entries in $\mathrm{Var}(\hat{\bm{\tau}}) $ ({\it i.e.}, sums of ratios of two higher-order polynomials of $z_{it}$). The objective of D-optimal design is the inverse of a $2(\ell+1)$-th order polynomial function of $z_{it}$. \cmtfinal{The objective functions of both A-optimal and D-optimal design are non-convex. Therefore, using the first order condition only is generally not sufficient to solve the global optimal solution.}

{\color{Black}

Optimal Solutions

We provide the optimality conditions for the T-optimal design. The optimality conditions disentangle the effect of different components in (ref) ({\it i.e.}, two-way fixed effects, observed and latent covariates) on the optimal treatment assignments. This problem is challenging in our setting with multiple units and periods for two reasons. First, different components can potentially affect the optimal design in both unit and time dimensions. Second, the effect of different components may be convoluted and interact with one another.

To build intuition, we start with the solution to a simple specification.

example[Two-way fixed effects and $\ell = 0$] Suppose $Y_{it} = \alpha_i + \beta_t + \tau_0 z_{it} + \varepsilon_{it}$, and and Assumptions (ref), (ref), and (ref) hold. Then the objective function in (ref) equals to \begin{align} \mathrm{Prec}(\hat{\tau}_0; A) = \frac{N}{\sigma_\varepsilon^2 } \left[ - 2 \bm{b}_T^\top \bm{\omega} -\bm{\omega}^\top \*P_{\bm{1}_{T}} \bm{\omega} \right], \end{align} which is a quadratic and concave function of $\bm{\omega} = [\omega_t]_{t \in [T]}$ with $\omega_t = N^{-1} \sum_{i=1}^N z_{it}$, $\*P_{\bm{1}_{T}} = \*I_{T} - \bm{1}_{T} \bm{1}_{T}^\top/T$, and $\*b_T = [b_{t}]$ with $b_{t} = (T+1-2t)/T$. By solving the first-order condition, we can show that any treatment design satisfying ${\omega}^{\ast}_t = -b_t$ for all $t$ is optimal and maximizes the precision. In the optimal solution, the treated fraction is linear in time. This result is conceptually similar to those in lawrie2015optimal,girling2016statistical,li2018optimal that show the optimal treated fraction is linear in time, under a similar specification, but with $\alpha_i$ to be random effects. More intuition for the linear treated fraction is provided in Section (ref).

Next consider a more general specification with $\ell > 0$, but without $\*X_i$ or $\*u_i$. We can show the objective function in (ref) is still quadratic and concave in $\bm{\omega}$, but takes a more complicated form (see Lemma (ref) in Section (ref)), leading to nonlinear optimal treated fraction in time. Specifically in Theorem (ref) we show that in the optimal solution, the unit average of $z_{it}$ satisfies

equation[equation omitted — 406 chars of source]

where $a^{(\ell)}$ is a vector of length $\ell - \lfloor \ell/2 \rfloor$ defined in Section (ref). $\omega_{\ell,t}^\ast$ has five stages and follows an $S$-shaped curve in time $t$: stage 1: all units are under control; stage 2: $\omega_t^\ast$ grows non-linearly in time (because of $a^{(\ell)}$); stage 3: $\omega_t^\ast$ grows linearly in time; stage 4: $\omega_t^\ast$ grows non-linearly in time again; stage 5, all units are under treatment. The optimal solution is symmetric with respect to the center ({\it i.e.}, $\left((T+1)/2, 0 \right)$). Figure (ref) demonstrates $\omega_{\ell,t}^\ast$ for various $\ell$ in a $T=12$ period problem. More examples of $\omega_{\ell,t}^\ast$ will be provided in Section (ref).

Furthermore, suppose the specification contains at least one of the $\*X_i$ or $\*u_i$ components. A commonly used variance reduction approach in cross-sectional studies is to balance covariates, so that treated and control units are comparable when the two groups have similar covariates ({\it e.g.}, imbens2015causal).\footnote{There is a strand of literature in operations research to use discrete optimization to achieve covariate balancing nikolaev2013balance,bertsimas2015power,bertsimas2019covariate,bhat2019near, and to use stratified sampling to increase power fox2000separability,mulvey1983multivariate.} Similar intuition carries over to the panel setting. We show in Theorem (ref) below that the optimal design balances covariates groups, where groups are defined as the set of units with the same initial treatment time.

We state the first of our main theorems using the following solution concepts: Let $\mathbb{A}_{\ell} = \big\{ A: N^{-1} \sum_{i = 1}^N \boldsymbol{1}_{A_i \leq t} = (1+\omega_{\ell,t}^\ast)/2,\, \forall t\big\}$ be the set of designs satisfying (ref). When $d_x > 0$, let $\mathbb{A}_{\*X} = \big\{ A: N^{-1} \sum_{i = 1}^N \*X_i \boldsymbol{1}_{A_i = t} = \bm{0}_{d_x},\, \forall t \in \{2,\cdots,T\}\big\}$ be the set of designs satisfying covariate balancing conditions (suppose rows in $\*X = [\*X_i]_{i \in [N]}$ are centered). When $d_x = 0$, let $\mathbb{A}_{\*X} = \mathbb{A}^N$ and no conditions need to be imposed related to $\*X_i$. $\mathbb{A}_{\*U}$ is defined similarly as $\mathbb{A}_{\*X}$. }

theorem[Optimality Conditions] Suppose Assumptions (ref), (ref), and (ref) hold, $\hat{\bm{\tau}} $ is \cmtrx{estimated from the infeasible GLS} with $\*W \propto \big(\*U \*\Sigma_v \*U^\top + \sigma_\varepsilon^2 \*I_N \big)^{-1} $, rows in $\*X$ and $\*U$ are centered $($i.e., $\sum_{i = 1}^N \*X_i = \bm{0}_{d_x}$ and $\sum_{i = 1}^N \*u_i = \bm{0}_{d_u}$$)$, $\*X$ and $\*U$ are orthogonal $($i.e., $\sum_{i = 1}^N \*X_i \*u_i^\top = \mathbf{0}_{d_x \times d_u}$$)$, $\*\Sigma_v = \sigma_\varepsilon^2 \cdot \*I_{d_u}$, and $T > (\ell^3+13\ell^2+7\ell+3)/(8\ell)$. $A$ is an optimal treatment design if $A \in \mathbb{A}_{\mathrm{opt}} $, where \begin{equation} \mathbb{A}_{\mathrm{opt}} = \mathbb{A}_{\ell} \cap \mathbb{A}_{\*X} \cap \mathbb{A}_{\*U}\, . \end{equation}
figure[figure omitted — 302 chars of source]

Theorem (ref) provides the optimality conditions for the T-optimal design when $\tau_0, \cdots, \tau_\ell$ are estimated from GLS using specification (ref). The presence of $\alpha_i$ and $\beta_t$ makes the optimal treated fraction grow gradually over time, and the growth rate is determined by $\ell$. We further illustrate how the optimal treated fractions depend on $\alpha_i$, $\beta_t$, and $\ell$ in Section (ref). The presence of $\*X_i$ and $\*u_i$ imposes additional covariate balancing conditions, that can be satisfied if units in each stratum satisfy treated fraction conditions. Here strata are defined as groups of units with the same covariate value. When covariates are discrete-valued, Example (ref) below provides a solution that satisfies the optimality conditions. For general cases, we provide guidance on choosing a design based on Theorem (ref) in Section (ref). \cmtrx{

example[Discrete $\*X_i$ and $\*u_i$] Suppose $(\*X_i, \*u_i)$ is discrete and can only take $G$ values for a finite $G > 0$, denoted as $\{(x_g, u_g)\}_{g=1}^G$, each with positive probability. Then any $A$ in $\mathbb{A}^{\mathrm{disc}}_{\mathrm{opt}}$ defined below is in $\mathbb{A}_{\mathrm{opt}}$ \begin{equation} \mathbb{A}^{\mathrm{disc}}_{\mathrm{opt}} =\bigg\lbrace A: \frac{1}{|\mathcal{O}_g|} \sum_{i \in \mathcal{O}_g} \boldsymbol{1}_{A_i \leq t} = \frac{1+\omega_{\ell,t}^\ast}{2} \,\,\,\, \forall t, g \bigg\rbrace\,, \end{equation} where $\mathcal{O}_g = \{i: (\*X_i, \*u_i) = (x_g, u_{g})\}$.

}

remark[Assumptions in Theorem (ref)] The assumptions in Theorem (ref) are non-restrictive for the following four reasons. First, the assumption that rows in $\*X$ and $\*U$ are demeaned can be satisfied by projecting the mean onto the unit fixed effects $\bm{\alpha}$. Second, the assumption that $\*X$ and $\*U$ are orthogonal can be satisfied by applying the Gram-Schmidt procedure ({\it i.e.}, QR decomposition) to $\begin{bmatrix} \*X & \*U \end{bmatrix}$ (which is possible as $\bm{\theta}$ is unknown and unrestricted). Third, $\*\Sigma_v = \sigma_\varepsilon^2 \cdot \*I_{d_u}$ is essentially an identification assumption, so that we can uniquely identify $\*u_i$ and $\*v_t$.\footnote{In other words, for arbitrary $\*u_i$ and $\*v_t$, we can right multiply $\*u_i^\top$ by $\sigma_\varepsilon \cdot \*\Sigma_v^{-1/2}$ and left multiply $\*v_t$ by $\sigma_\varepsilon \cdot \*\Sigma_v^{-1/2}$ so that $v_t$ has variance $\sigma_\varepsilon^2 \cdot \*I_{d_u}$ (conditions $\sum_{i = 1}^N \*u_i = \bm{0}_{d_x}$ and $\sum_{i = 1}^N \*X_i \*u_i^\top = \mathbf{0}_{d_x,d_u}$ stay valid after this manipulation).} Fourth, the fundamental structure of our problem does not change with these assumptions because $\mathrm{Prec}(\hat{\bm{\tau}})$ is a quadratic function of $z_{it}$ regardless of whether these assumptions are imposed or not.
remark[Assumption on $T$] The assumption $T > (\ell^3+13\ell^2+7\ell+3)/(8 \ell) $ is a sufficient condition to show $\omega_{\ell,t}^\ast$ is monotonic, but is not a necessary condition and can potentially be relaxed, based on the numerical solution of treated fractions when this assumption is violated.
remark[Magnitude of treatment effects] Note that $\mathbb{A}_{\mathrm{opt}} $ and $\mathbb{A}^{\mathrm{disc}}_{\mathrm{opt}}$ do not depend on the value of $\bm{\tau}$. This is because the estimation error $\hat{\bm{\tau}} - \bm{\tau}$ from GLS does not depend on the value of $\bm{\tau}$ under the specification (ref). Therefore, $\mathrm{tr} \big(\mathrm{Prec}(\hat{\bm{\tau}}) \big)$ does not depend on the value of $\bm{\tau}$. See (ref) for an example.

Optimal Treated Fractions

figure[figure omitted — 1,040 chars of source]

{\color{Black}

To provide more intuition on the optimal treated fractions, we start with the specification with $\ell = 0$, and with either $\alpha_i$ or $\beta_t$, but not both.

example[Time Fixed Effects Only and $\ell = 0$] Suppose $Y_{it} = \beta_t + \tau_0 z_{it} + \varepsilon_{it}$ and Assumptions (ref), (ref), and (ref) hold. In this case, $\mathrm{Prec}(\hat{\tau}_0) = {N}/{\sigma_\varepsilon^2 } \cdot \big[T - \sum_{t=1}^T \omega_t^2 \big]$. Any treatment design is optimal if it satisfies $\omega^\ast_t = 0$ for all $t$ and $z_{it} \leq z_{i,t+1}$ for all $i$ and $t$, that is, assigns treatment to 50% units, while the others functioning as the control, for all time periods.
example[Unit Fixed Effects Only and $\ell = 0$] Suppose $Y_{it} = \alpha_i + \tau_0 z_{it} + \varepsilon_{it}$ and Assumptions (ref), (ref), and (ref) hold. In this case, $\mathrm{Prec}(\hat{\tau}_0) = {T}/{\sigma_\varepsilon^2 } \cdot \big[N - \sum_{i=1}^N \zeta_i^2 \big]$, where $\zeta_i = T^{-1} \sum_{t=1}^T z_{it}$. Any treatment design is optimal if it satisfies $\omega^\ast_t = -1$ for all $t < (T+1)/2$ and $\omega^\ast_t = 1$ for all $t > (T+1)/2$, that is, allocates treatment to all units at halftime.\footnote{For odd $T$, units can be either treated or untreated at $t = (T+1)/2$.}

}

If the specification has both $\alpha_i$ and $\beta_t$, as in Example (ref), the optimal treated fractions are intuitively in between those in Examples (ref) and (ref). See Figure (ref) for the visualization of optimal treated fractions under various specifications.

The same intuition carries over to $\ell > 0$ and the precision matrix $\mathrm{Prec}(\hat{\bm{\tau}})$ is quadratic and concave in $\omega_t$. The two examples below provide the analytic expression of $\omega^\ast_{\ell,t}$ for $\ell = 1$ and $2$. We provide the expression of $\omega_{\ell,t}^\ast$ for $\ell = 3$ in Example (ref) of Section (ref).

example[$\ell = 1$] In Theorem (ref), $\omega^\ast_{\ell,t}$ equals $-1 + {2(t-1)}/{(T-1)}$ for all $t$.
example[$\ell = 2$] In Theorem (ref), $\omega^\ast_{\ell,t}$ is determined by, $\omega_{\ell,1}^\ast = -1$, $\omega_{\ell,2}^\ast = -1 + {2}/{(2T-5)}$, $\omega_{\ell,t}^\ast = -1 + {(2t-3)}/{(T-2)}$ for $t = 3, \cdots, T-2$, $\omega_{\ell,T-1}^\ast = 1 - {2}/{(2T-5)}$, and $\omega_{{\ell,T}}^\ast = 1$.

{\color{Black} \paragraph{Intuition for the $S$-shaped curve of $\omega_{\ell,t}^\ast$.} The objective function in (ref) is a sum of $\ell+1$ quadratic functions, one function for each $\mathrm{Prec}(\hat{\tau}_j)$ with $j \in \{0,1,\cdots,\ell\}$. Let $\omega_{\ell,t}^\ast(j)$ be the unit average of $z_{it}$ in the solution that maximizes $\mathrm{Prec}(\hat{\tau}_j)$ at time $t$, when the duration of lagged effects is $\ell$. This means $\omega^\ast_{\ell,t}$, defined in (ref), can be written as a convex combination of $\omega_{\ell,t}^\ast(j)$ for all $j$ and $t$. Intuitively, since each $\omega_{\ell,t}^\ast(j)$ is chosen to maximize the precision of a single parameter $\hat{\tau}_j$, it should grow linearly in $t$, like the $\ell=0$ case, except during the first $\ell$ or the last $\ell$ periods that it is truncated to $0$ or $1$ respectively (see Figure (ref)). Equipped with this observation, it is easy to see that $\omega^\ast_{\ell,t}$ grows linearly in $t$ in the middle (the third stage) where all $\omega^\ast_{\ell,t}$ grow linearly. $\omega^\ast_{\ell,t}$ is constant in stage one and five, where all $\omega^\ast_{\ell,t}$ are constant. Moreover, $\omega^\ast_{\ell,t}$ has a non-linear growth in stage two and four, because some of $\omega_{\ell,t}^\ast(j)$ are truncated to $0$ or $1$. }

figure[figure omitted — 747 chars of source]

Choosing A Treatment Design

{\color{Black} Based on the optimality conditions in Theorem (ref), we provide Algorithm (ref) in Section (ref) to choose a treatment design.} In this subsection, we outline some general guidance on choosing a treatment design. When the specification does not have covariates ($d_x = d_u = 0$) and $\mathbb{A}_{\ell} $ is not empty, we can randomly choose one from $\mathbb{A}_\ell$ with equal probability. When the specification has observed, discrete-valued covariates ($d_x \neq 0$) and $\mathbb{A}^{\mathrm{disc}}_{\mathrm{opt}}$ is not empty, we can randomly choose one from $\mathbb{A}^{\mathrm{disc}}_\ell$ with equal probability. The random sampling can ensure that the treatment design balances the relevant covariates for the outcomes that are omitted in the specification (ref) and design of experiments hayes2017cluster.

Note that $\mathbb{A}^{\mathrm{disc}}_{\mathrm{opt}}$ can be empty, if $|\mathcal{O}_g|(1+\omega_{\ell,t}^\ast)/(2T)$ is not an integer for each stratum (similarly for $\mathbb{A}_{\ell} $). For this case, Algorithm (ref) rounds $|\mathcal{O}_g|(1+\omega_{\ell,t}^\ast)/(2T)$ to the nearest integer to obtain a feasible design. As shown in Proposition (ref) in Section (ref), the value of $\mathrm{tr} \big(\mathrm{Prec}(\hat{\bm{\tau}}) \big)$ evaluated at this feasible solution is within a factor of $1 + O \left( {1}/{N^2} \right)$ of the optimal value of $\mathrm{tr} \big(\mathrm{Prec}(\hat{\bm{\tau}}) \big)$.

{\color{Black} If observed covariates are continuous, then we can partition units into a small number of strata based on the observed covariate values, and then randomly choose a design that satisfies the treated fraction conditions for each stratum (possibly with rounding). Prior work has suggested keeping the number of strata small to avoid over-stratification kernan1999stratified, because over-stratification may lower the precision of the estimated treatment effects de2008consequences. }

{\color{Black} If there are latent covariates, we suggest using the historical control data for the same set of units for the design of experiments. We can improve precision by using historical data to estimate $\*u_i$, partitioning units into strata ({\it e.g.}, by spectral clustering), and choosing a treatment design based on the estimated $\*u_i$. See Section (ref) for more details. }

remark[Conditional Average Treatment Effect (CATE)] Suppose we are interested in using observable sources of heterogeneity $\*X_i$ to assess the heterogeneity in treatment effects (robinson1988root,wager2018estimation among others). In this case, we may be interested in conditional average instantaneous and lagged effects, denoted by $\tau_j(\*X_i)$. If $\*X_i$ is discrete, the designs satisfying (ref) can maximize the precision of $\hat{\tau}_j(\*X_i)$ for all $j$ and $\*X_i$, following that the optimality conditions for $\omega_{\ell,t}^\ast$ are satisfied within each stratum.

{\color{Black}

remark[Discussion on Adaptive Design] The treatment assignments studied in this section are chosen and fixed before the experiment starts. For non-adaptive experiments, we do not pursue adaptive designs, where treatment assignments for subsequent periods can vary during the experiment, for the following reason. From Theorem (ref), the information that matters for the treatment assignments but is unknown before the experiment starts is $\*u_i$. However, the estimation of $\*u_i$ at the beginning of the experiment can be quite noisy, and we do not expect noisy estimates of $\*u_i$ can substantially help the design. See Section (ref) for further elaboration on this.

}

\setcounter{equation}{0} {\thesection.\arabic{equation}} \oldsection{Adaptive Experiments}

Consider an experiment that is allowed to run $T_{\max}$ periods. The experimenter seeks to achieve a certain precision objective by the end of the experiment. After observing some data during the experiment, the experimenter may find that achieving this objective does not need $T_{\max}$ periods' of experimental data. In this case, the experimenter may want to terminate the experiment early in order to reduce the cost of the experiment.

In this section, we study the design and analysis of adaptive experiments, where $N$ is fixed, but the experiment duration can vary due to the early termination. Let $\tilde{T} \in [T_{\max}]$ be the duration of the adaptive experiment, which is a random variable and unknown before the adaptive experiment starts. At any time period $t$, the experimenter collects data and decides whether to terminate the experiment. If so, the experiment stops, and the realization of $\tilde{T}$ is $t$; otherwise, the experimenter makes treatment decisions for time $t + 1$. For ease of understanding, we present our algorithm and results based on the following simple specification of the observed outcome of unit $i$ at time $s$

equation[equation omitted — 166 chars of source]

For simplicity, we denote instantaneous effect as $\tau$ instead of $\tau_0$ in this section. Later, we discuss how our algorithm and results can be generalized to the specification with $\ell > 0$ and with $\*X_i$ and $\*u_i$.

Motivated by the objective of maximizing precision in (ref), we consider the following criterion to terminate the experiment if the precision exceeds a target threshold $c$ at time $t$

equation[equation omitted — 296 chars of source]

\cmtrx{where the expression of $\mathrm{Prec}(\hat{\tau}; \bm{\omega})$ follows from Example (ref) in the case of $\hat{\tau}$ estimated from a panel with $N$ units and $t$ periods. In this section, we parametrize the precision by $\bm{\omega}$ and we optimize over $\bm{\omega}$, as the precision only depends on $\bm{\omega}$ but not whom to treat under specification (ref).} This precision-based rule is equivalent to the variance-based rule to terminate the experiment when $\mathrm{Var}(\hat{\tau}) \leq 1/c$. This type of termination rules has been used by others in different settings, such as in simulations and in sequential testing on sequentially arrived units without time effects (chow1965asymptotic,glynn1992asymptotic,singham2012finite among others).

There are three key technical challenges in designing and analyzing the experiments that can be terminated early. The first challenge concerns adaptively \cmtfinal{choosing the fraction of treated units per period}. Recall from Example (ref), that $\omega_s^\ast =(2s-1-\tilde{T})/\tilde{T}$. But, since $\tilde{T}$ is unknown before the experiment starts, or even during the experiment, \cmtfinal{choosing the optimal fraction of treated units} is non-trivial. To address this challenge, we aim to adaptively improve the treatment decisions as we gather more information about $\tilde{T}$ during the experiment.

The second challenge concerns implementing the termination rule. As long as we can estimate the critical unknown parameter $\sigma_\varepsilon^2$ in (ref), we can determine whether to stop the experiment. There are two main difficulties in this task. The first is to have a valid implementation of the termination rule on adaptively collected data. Here, valid implementation means that the precision indeed exceeds threshold $c$ when the experiment terminates. The second is to do it in a way that the next challenge (obtaining valid post-experiment inference of $\tau$) can be manageable. Early stopping complicates post-experiment inference because the same data is used to make the decision about stopping and to estimate treatment effects, leading to the well-known bias that can arise when adaptive tests determine whether to terminate experiments (johari2017peeking among others).

The third challenge concerns efficient estimation and inference for $\tau$, post-experiment. Based on the adaptive nature of collected data and the implementation of experiment termination rule, we seek to choose a consistent and efficient estimator for $\tau$ that uses as many observations as possible.

We propose the Precision-Guided Adaptive Experiment (PGAE) algorithm in Section (ref) to simultaneously tackle these three challenges. PGAE combines ideas from dynamic programming and sample splitting. In Section (ref), we prove statistical consistency and asymptotic normality of $\hat\tau$ and $\widehat{\sigma^2}$, estimated by PGAE, paving the way for valid statistical inference for $\tau$.

Estimators

To start, we first review two existing estimators and then propose a new estimator. All three estimators are extensively used in PGAE. Suppose these estimators use the data of units in a set $\mathcal{S}$ over $t$ periods collected so far, where $t$ is small, but set size $|\mathcal{S}|$ can be large. In this subsection, we sub-index the estimators by $\mathcal{S}$ and $t$ to refer to the data used in the estimators.

The first is the within estimator for $\tau$ wallace1969use. The within estimator of $\tau$, denoted by $\hat{\tau}_{\mathcal{S},t}$, regresses $\dot{Y}_{is}$ on $\dot{z}_{is}$ based on the specification $\dot{Y}_{is} = \tau \dot{z}_{is} + \dot{\varepsilon}_{is}$, where for any variables $\{x_{is}\}_{(i,s) \in \mathcal{S} \times [t]}$ ({\it e.g.}, $Y_{is}$ and $z_{is}$), the notation $\dot{x}_{is}$ denotes the within transformed $x_{is}$ and is defined as

equation[equation omitted — 122 chars of source]

in which $\bar{x}_{i \cdot}$, $\bar{x}_{\cdot s}$, and $ \bar{x}$ are averages of $x_{is}$'s over time periods, units, and both of them, respectively. The within estimator is an efficient estimation approach that does not need to estimate $\alpha_i$ and $\beta_s$, but produces the same estimate of $\tau$ as OLS based on (ref) that estimates $\alpha_i$ and $\beta_s$.\footnote{Regressing $Y_{is}$ on $z_{is}$ and unit and time dummies is the same as GLS with weight matrix $\*W \propto \*I_N$ in (ref). The within estimator is also called Least-Squares Dummy Variable (LSDV) estimator. } As shown in Lemma (ref), $\hat{\tau}_{\mathcal{S},t}$ is consistent and asymptotically normal for any finite $t$, when the set size $|\mathcal{S}|$ is large.

The second is the plug-in estimator for $\sigma^2_\varepsilon$, which is used in experiment termination and post-experiment inference, and takes the form of

equation[equation omitted — 251 chars of source]

The factor $1/(t-1)$ is for finite $t$ correction. As shown in Lemma (ref), $\widehat{\sigma^2}_{\mathcal{S},t}$ is consistent and asymptotically normal for any finite $t$.

The third is a new estimator for the variance of $\varepsilon_{is}^2$, that is, $\xi^2_\varepsilon \coloneqq \+E[(\varepsilon^2_{is} - \sigma^2_\varepsilon)^2]$, which is used to quantify the uncertainty in our estimator for $\sigma^2_\varepsilon$, and takes the form of

align[align omitted — 562 chars of source]

In this estimator, both the correction multiplier and correction term are used to de-bias the plug-in estimator when $t$ is finite. The plug-in estimator is biased, because $(\dot{y}_{is} - \hat{\tau}_{\mathcal{S},t} \cdot \dot{z}_{is})^2$ is not an unbiased estimator of $\sigma^2_{is}$ for each $i$ and $s$, and the bias is squared in the estimation of $\xi_\varepsilon^2$, which cannot be averaged out over $i$. Figure (ref) in Section (ref) visualizes the bias of the plug-in estimator and shows that $\widehat{\xi^2}_{\mathcal{S},t}$ can correct for the bias in finite samples. Lemma (ref) shows that $\widehat{\xi^2}_{\mathcal{S},t}$ is consistent for any finite $t$.

Precision-Guided Adaptive Experiment (PGAE) via Dynamic Programming

PGAE simultaneously addresses the three challenges introduced at the beginning of Section (ref), which is feasible through a careful partitioning of units into disjoint subsets, with each set functioning for a different purpose. Specifically, we partition units into three mutually disjoint sets $\mathcal{S}_{\mathrm{ntu}}$, $\mathcal{S}_{\mathrm{atu},1}$ and $\mathcal{S}_{\mathrm{atu},2}$ that represent a set of non-adaptive treatment units (NTU) and two sets of adaptive treatment units (ATU), respectively. See Figure (ref) for an illustration. Let $p_{\mathrm{ntu}} = |\mathcal{S}_{\mathrm{ntu}}|/N$, $p_{\mathrm{atu},1} = |\mathcal{S}_{\mathrm{atu},1}|/N$, and $p_{\mathrm{atu},2} = |\mathcal{S}_{\mathrm{atu},2}|/N$ be the fractions of units in these three sets. Clearly, $p_{\mathrm{ntu}} + p_{\mathrm{atu},1} + p_{\mathrm{atu},2} = 1$. The sets are selected such that $p_{\mathrm{ntu}}$ is small and $p_{\mathrm{atu},1} = p_{\mathrm{atu},2}$.

Before the experiment starts, we initialize the treatment designs of $\mathcal{S}_{\mathrm{ntu}}$, $\mathcal{S}_{\mathrm{atu},1}$ and $\mathcal{S}_{\mathrm{atu},2}$ by the optimal design of a $T_{\max}$-period non-adaptive experiment, which is the solution without early stopping. Then the average of $z_{is}$ over $i$ in each set satisfies $\omega_{\mathrm{bm},s} = (2s - 1 - T_{\max})/{T_{\max}}$ for all $s\in[T_{\max}]$ (with rounding if necessary), per Example (ref). The initial design also serves as the “benchmark” design for the adaptive experiment. The treatment design of NTU does not change during the experiment, stays equal to $\bm{\omega}_{\mathrm{bm}}$, and observed data from NTU is used to update the treatment assignments of ATU for subsequent periods, specifically, to improve upon $\bm{\omega}_{\mathrm{bm}}$. One of the two ATU sets ($\mathcal{S}_{\mathrm{atu},1}$) is used to estimate $\sigma_\varepsilon^2$, and decide whether to terminate the experiment. The other ATU set ($\mathcal{S}_{\mathrm{atu},2}$) provides another estimate of $\sigma_\varepsilon^2$ for post-experiment inference. As we will see in Section (ref), partitioning units into three sets is crucial for decoupling possible correlations due to data reuse and hence obtaining valid statistical inference.

figure[figure omitted — 162 chars of source]

Next, we illustrate how PGAE addresses each of the three challenges.

\paragraph{1. Choosing a treatment design.} {\color{Black} In order to update treatment decisions for ATU, we need to construct a belief about experiment stopping time $\tilde{T}$, using NTU. If we would have known $\sigma_\varepsilon^2$, then we would exactly know the minimum stopping time $\tilde{T}$ such that the precision (using $\bm{\omega}_{\mathrm{bm}}$) is bigger than $c$ in (ref). However, $\sigma_\varepsilon^2$ is unknown in practice, but it is possible to construct \cmtfinal{the belief} distribution of $\sigma_\varepsilon^2$, which will be explained in detail below. Then we can draw samples of $\sigma^2$ from the \cmtfinal{belief} distribution of $\sigma_\varepsilon^2$. For each sampled $\sigma^2$, we plug it into the precision expression and find the minimum duration that satisfies the stopping rule (ref).\footnote{Formally, we find $T$ using the following approach. Let $\sigma_0^2$ be a sampled $\sigma^2$ and let $T_0 = \min \left\{t: Nt / \sigma_0^2 \cdot g_\tau(\bm{\omega}, t) \geq c \right\}$ be the minimum duration such that the precision exceeds threshold $c$. If $T_0 \in \{t+1, t+2, \cdots, T_{\max}\}$, then we set $T$ as $ T_0$; otherwise, if $T_0 > T_{\max}$, then we set $T$ as $T_{\max}$; otherwise, if $T_0 \leq t$, then we set $T$ as $t+1$. } The minimum duration may vary with the sampled $\sigma^2$. By repeatedly sampling $\sigma^2$ and finding the minimum duration, we can obtain an empirical distribution of $\tilde{T}$, denoted as $P_t(\tilde{T})$. See the pseudocode of the helper function estimate_belief() of PGAE in Section (ref) for more details.

Here is our approach to constructing a \cmtfinal{belief} distribution of $\sigma_\varepsilon^2$. We first estimate $\sigma_\varepsilon^2$ from formula (ref) using data of NTU for $t$ periods. Let $\widehat{\sigma^2}_{\mathrm{ntu},t}$ be the estimator. Next, we quantify the uncertainty in $\widehat{\sigma^2}_{\mathrm{ntu},t}$. From Lemma (ref) below, $\widehat{\sigma^2}_{\mathrm{ntu},t}$ is consistent and asymptotically normal with

align[align omitted — 249 chars of source]

where $\hat{\xi}_{\mathrm{ntu},t}^{\dagger} = \big[\widehat{\xi^2}_{\mathrm{ntu},t} + (\widehat{\sigma^2}_{\mathrm{ntu},t})^2/(t-1) \big]^{1/2}$ and $\widehat{\xi^2}_{\mathrm{ntu},t}$ is estimated from (ref) using data of NTU. Let the left-hand side of (ref) be $w$. We have $\sigma_\varepsilon^2 = \widehat{\sigma^2}_{\mathrm{ntu},t} - w \cdot {\hat{\xi}^\dagger_{\mathrm{ntu},t}}/{\sqrt{|\mathcal{S}_\mathrm{ntu}| \cdot t}}$, and we repeatedly use this formula with samples of $w$ from $\mathcal{N}(0,1)$ to obtain a \cmtfinal{belief} distribution of $\sigma_\varepsilon^2$.

We use a dynamic program to solve the treatment decisions for ATU. In this dynamic program, the state variable is the treated fractions up to time $t$, i.e., $\bm{\omega}_{\mathrm{atu},1:t}$, and the decision variable is the treated fraction for the next period $\omega_{t+1}$.\footnote{{\color{Black} Example 4.3.4 of bertsekas2012dynamic also uses a dynamic program to select a threshold for terminating a sequential hypothesis testing problem, but we use a dynamic program to find the optimal design in subsequent time periods.}} The payoff in intermediate time periods is zero, and conditional on the realization of $\tilde{T}$, the terminal cost is $- \tilde{T} \cdot g_\tau\left((\bm{\omega}_{\mathrm{atu},1:t}, \bm{\omega}_{(t+1):\tilde{T}}), \tilde{T}\right)$, which is equivalent to our objective of maximizing precision. Here, in the dynamic program, we aim to solve $\omega_{t+1}$ that minimizes the expected terminal cost. The expectation is taken with respect to the random experiment duration $\tilde{T}$, and we use the empirical distribution $P_t(\tilde{T})$ learned at time $t$ when taking the expectation. Specifically, we solve $\omega_t$ by minimizing \[ \+E_{\tilde{T} \sim P_t(\tilde{T})} \left[ - \tilde{T} \cdot g_\tau\left((\bm{\omega}_{\mathrm{atu},1:t}, \bm{\omega}_{(t+1):\tilde{T}}), \tilde{T}\right) \right]\,, \] subject to the constraint $\omega_{\mathrm{atu},t} \leq \omega_{t+1} \leq \omega_{t+2} \leq \cdots \leq \omega_{T_{\max}} \leq 1$.\footnote{Only $\omega_{t+1}$ is the decision variable. $\omega_{t+2}, \cdots, \omega_{T_{\max}} $ are not decision variables, whose values can vary with $\tilde{T}$.} Section (ref) provides more discussion on this dynamic program, and more details about the dynamic program can be found in the pseudocode of the helper function update_treatment_design() of PGAE in Section (ref).

}

\paragraph{2. Implementing the termination rule.} We use the data in $\mathcal{S}_{\mathrm{atu},1}$ to estimate $\sigma^2_\varepsilon$ via formula (ref). Let the estimator at time $t$ be $\widehat{\sigma^2}_{\mathrm{atu},1,t}$.\footnote{For notation simplicity, we denote the estimator as $\widehat{\sigma^2}_{\mathrm{atu},1,t}$, as opposed to $\widehat{\sigma^2}_{\mathcal{S}_{\mathrm{atu},1},t}$.} We then plug $\widehat{\sigma^2}_{\mathrm{atu},1,t}$ and $\bm{\omega}_{\mathrm{bm},1:t}$ into the termination rule (ref) to estimate precision and decide whether to terminate the experiment. Note that using $\bm{\omega}_{\mathrm{bm},1:t}$ tends to under-estimate the precision (see Proposition (ref) in Section (ref)), so that the experiment termination tends to be conservative. If our experiment terminates, we show in Proposition (ref) in Section (ref) that the true precision exceeds $c$ with a high probability.\footnote{Admittedly, it seems natural to use $\bm{\omega}_{\mathrm{atu},1:t}$ that may yield a more precise estimation of the precision. We do not pursue this route as it is harder to show that the post-experiment inference of $\tau$ is valid and the corresponding precision is indeed larger than $c$, given the complexity of the current proof.}

\paragraph{3. Efficient estimation and valid inference for $\tau$.} In PGAE, we estimate $\tau$ using all the data of NTU and ATU over $\tilde{T}$ periods. Let the estimator be $\hat{\tau}_{\mathrm{all},\tilde T}$. We show in Theorem (ref) below that $\hat{\tau}_{\mathrm{all},\tilde T}$ is efficient and achieves the optimal convergence rate. But we only use $\mathcal{S}_{\mathrm{atu},2}$ to estimate $\sigma^2_\varepsilon$ via formula (ref). Let the estimator be $\widehat{\sigma^2}_{\mathrm{atu},2,\tilde T}$, which is consistent, and we use $\widehat{\sigma^2}_{\mathrm{atu},2,\tilde T}$ to estimate $\mathrm{Var}(\hat{\tau})$. Note that the consistency of $\sigma^2_\varepsilon$ is sufficient for constructing valid confidence intervals for $\tau$. Therefore the efficiency (determined by the sample size) in the estimation of $\sigma^2_\varepsilon$ is less important. The only tradeoff is that the endpoints of confidence intervals are estimated less precisely, with a larger second-order error term.\footnote{The first-order error term comes from the estimation error of $\hat{\tau}_{\mathrm{all},\tilde T}$} This explains why we use only $\mathcal{S}_{\mathrm{atu},2}$ to estimate $\sigma^2_\varepsilon$.

In summary, PGAE takes advantage of all the data collected so far, and then jointly optimizes treatment assignments alongside the choice of whether to continue the experiment. The pseudocode for PGAE is shown in Algorithm (ref), with the pseudocode for the helper functions collected in Section (ref).

algorithm[algorithm omitted — 2,149 chars of source]

Analysis of the Algorithm

In this subsection, we present the asymptotic results of the estimated $\tau$ and $\sigma^2_\varepsilon$ in PGAE. These results serve two purposes. The first is to justify our approach in Section (ref) to construct a belief about $\sigma^2_\varepsilon$, and specifically to justify (ref). The second is to show that the outputs of PGAE ({\it i.e.}, $\hat{\tau}_{\mathrm{all},\tilde{T}}$ and $\widehat{\sigma^2}_{\mathrm{atu},2,\tilde{T}}$) can be used for valid statistical inference and hypothesis test of treatment effect.

{\color{Black} To start, we characterize the asymptotic properties of $\hat{\tau}_{\mathrm{ntu},t}$, $\widehat{\sigma^2}_{\mathrm{ntu},t}$ and $\widehat{\xi^2}_{\mathrm{ntu},t}$, when they are estimated from non-adaptive experimental data, in Lemma (ref). This lemma provides theoretical support for constructing the \cmtfinal{belief} distribution of $\sigma_\varepsilon^2$, and serves as a crucial intermediate step in characterizing asymptotic distributions for the outputs of PGAE.

lemmaSuppose Assumptions (ref) and (ref) hold. Under the specification (ref), suppose $\varepsilon_{is}$ is i.i.d. for any $i$ and $s$ with $\+E[\varepsilon_{is}] = 0$, $\+E[\varepsilon^2_{is}] = \sigma_\varepsilon^2$, $\+E[\varepsilon^3_{is}] = 0$, and $\+E[(\varepsilon^2_{is} - \sigma^2_\varepsilon)^2] = \xi^2_{\varepsilon}$. $\hat{\tau}_{\mathrm{ntu},t}$ and $\widehat{\sigma^2}_{\mathrm{ntu},t}$ are consistent. As $|\mathcal{S}_{\mathrm{ntu}}| \rightarrow \infty$, for any finite $t$, conditional on $Z_\mathrm{ntu}$, we have \[\sqrt{|\mathcal{S}_{\mathrm{ntu}}|} \left(\begin{bmatrix} \hat{\tau}_{\mathrm{ntu},t} \\ \widehat{\sigma^2}_{\mathrm{ntu},t} \end{bmatrix} - \begin{bmatrix} \tau \\ \sigma^2_\varepsilon \end{bmatrix}\right) \stackrel{d}{\longrightarrow } \mathcal{N} \left(\begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} \sigma_\varepsilon^2/(t \cdot g_\tau(\bm{\omega}_{\mathrm{ntu}, 1:t}, t)) & 0 \\ 0 & \xi_{\varepsilon,t}^{\dagger2}/t \end{bmatrix} \right), \] where $\xi_{\varepsilon,t}^{\dagger2} = \xi^2_\varepsilon + 2 \big(\sigma_\varepsilon^2\big)^2/(t-1)$. Furthermore, $\sqrt{|\mathcal{S}_{\mathrm{ntu}}|} \big(\widehat{\xi^2}_{\mathrm{ntu},t} - \xi_\varepsilon^2\big) = O_p(1)$.

Since both $\widehat{\sigma^2}_{\mathrm{ntu},t}$ and $\widehat{\xi^2}_{\mathrm{ntu},t}$ are consistent, $\widehat{\xi^{\dagger 2}}_{\mathrm{ntu},t} = \widehat{\xi^2}_{\mathrm{ntu},t} + 2 \big(\widehat{\sigma^2}_{\mathrm{ntu},t} \big)^2/(t-1)$ is a consistent estimator of $\xi_\varepsilon^{\dagger2}$. From Slutsky's theorem, the asymptotic distribution in (ref) holds.

Interestingly, $\hat{\tau}_{\mathrm{ntu},t}$ and $\widehat{\sigma^2}_{\mathrm{ntu},t}$ are asymptotically independent, even though we use $\hat{\tau}_{\mathrm{ntu},t}$ in the estimation of $\widehat{\sigma^2}_{\mathrm{ntu},t}$. The reason of asymptotic independence is as follows. The estimation error of $\hat{\tau}_{\mathrm{ntu},t}$ is a weighted average of $\varepsilon_{is}$ over $i$ and $s$. The leading term in the estimation error of $\widehat{\sigma^2}_{\mathrm{ntu},t} $ is the sum of a weighted average of $\varepsilon^2_{ju} - \sigma_\varepsilon^2$ over $j$ and $u$ and a weighted average of $\varepsilon_{ju} \varepsilon_{jv} $ over $j,u,v$ with $u \neq v$. Note that $\varepsilon_{is}$ is uncorrelated with both $\varepsilon^2_{ju} - \sigma_\varepsilon^2$ and $\varepsilon_{ju} \varepsilon_{jv} $ for all $i,j,s,u,v$, because $\varepsilon_{is}$ is i.i.d. in $i$ and $s$ and has zero first and third moments. Therefore, the leading terms in the estimation errors of $\hat{\tau}_{\mathrm{ntu},t}$ and $\widehat{\sigma^2}_{\mathrm{ntu},t} $ are uncorrelated. For the non-leading terms, they are at a small order of magnitude and do not contribute to the asymptotic covariance. As $\hat{\tau}_{\mathrm{ntu},t}$ and $\widehat{\sigma^2}_{\mathrm{ntu},t}$ are jointly asymptotically normal, the zero asymptotic correlation implies asymptotic independence. In Section (ref), we demonstrate the finite sample properties of Lemma (ref) and asymptotic independence between $\hat{\tau}_{\mathrm{ntu},t}$ and $\widehat{\sigma^2}_{\mathrm{ntu},t}$.

remarkLemma (ref) holds for any finite $t$. In fact, when $t$ grows to infinity, the problem is simpler, because we can show $\mathrm{plim}_{t\rightarrow \infty} \xi^{\dagger 2}_{\varepsilon,t} = \xi^2_\varepsilon$, and the plug-in estimator of $\xi^2_\varepsilon$ mentioned in formula (ref) is consistent. In this section, we focus on the challenging case with a finite $t$, because we want to apply Lemma (ref) to the estimates on NTU early in the experiment ({\it i.e.}, $t$ is small).

}

Next we show that $\hat{\tau}_{\mathrm{all},\tilde{T}}$ and $\widehat{\sigma^2}_{\mathrm{atu},2,\tilde{T}}$ from PGAE can be used for valid post-experiment statistical inference and hypothesis testing for $\tau$. To show this, there are two critical steps: (a) show the asymptotic distribution of $\hat{\tau}_{\mathrm{all},\tilde{T}}$; (b) show that the asymptotic variance of $\hat{\tau}_{\mathrm{all},\tilde{T}}$ can be consistently estimated using $\widehat{\sigma^2}_{\mathrm{atu},2,\tilde{T}}$.

theoremSuppose Assumptions (ref) and (ref) hold, $p_{\mathrm{ntu}} \in (0,1)$ is a fixed number as $N$ grows. Under the specification (ref), suppose the assumptions about $\varepsilon_{is}$ in Lemma (ref) hold and $\varepsilon_{is}$ is bounded with a symmetric distribution around $0$. $\hat{\tau}_{\mathrm{all},\tilde{T}}$ and $\widehat{\sigma^2}_{\mathrm{atu},2,\tilde{T}}$ are consistent. As $N \rightarrow \infty$, \begin{align} \sqrt{N} \cdot \begin{bmatrix} \big(\tilde{T}g_\tau(\bm{\omega}_{\mathrm{all},1:\tilde{T}},\tilde{T})/\sigma_\varepsilon^2\big)^{1/2} \cdot \left( \hat{\tau}_{\mathrm{all},\tilde{T}} - \tau\right) \\ \big(\tilde{T} p_{\mathrm{atu},2}/\xi^{\dagger 2}_{\varepsilon,\tilde{T}} \big)^{1/2} \cdot \big(\widehat{\sigma^2}_{\mathrm{atu},2,\tilde{T}} - \sigma_\varepsilon^2 \big) \end{bmatrix} \stackrel{d}{\longrightarrow } \mathcal{N} \left(\bm{0}, I_2 \right) \, . \end{align}

{

\cmtrx{$\hat{\tau}_{\mathrm{all},\tilde{T}}$ is consistent for $\tau$ with the optimal convergence rate $\sqrt{N}$. This result is not obvious because $\mathcal{S}_{\mathrm{ntu}}$, $\mathcal{S}_{\mathrm{atu},1}$ and $\mathcal{S}_{\mathrm{atu},2}$ are all used to estimate $\tau$, which could potentially lead to two sources of bias in the estimation of $\tau$. The first is that $\hat{\tau}_{\mathrm{all},\tilde{T}}$ depends on adaptive treatment designs, and the choice of adaptive designs depend on $\widehat{\sigma^2}_{\mathrm{ntu},t}$ and $\widehat{\xi^2}_{\mathrm{ntu},t}$, and therefore on $\varepsilon_{is}$ for $i \in \mathcal{S}_{\mathrm{ntu}}$. The second is that $\hat{\tau}_{\mathrm{all},\tilde{T}}$ depends on termination time $\tilde{T}$, where $\tilde{T}$ depends on $\widehat{\sigma^2}_{\mathrm{atu},1,t}$ and therefore on $\varepsilon_{is}$ for $i \in \mathcal{S}_{\mathrm{atu},1}$. Both sources may lead to the violation of the commonly made exogeneity assumption to show consistency ({\it i.e.}, asymptotic conditional mean of $\varepsilon_{it}$ is zero). With a careful analysis in Lemma (ref), we show this is not the case, and the asymptotic conditional mean is still zero. This is mainly because $\widehat{\sigma^2}_{\mathrm{ntu},t}$, $\widehat{\xi^2}_{\mathrm{ntu},t}$ and $\widehat{\sigma^2}_{\mathrm{atu},1,t}$ are all even moments of $\varepsilon_{is}$. With a symmetric distribution of $\varepsilon_{is}$ around $0$, conditioning on the even moments does not change the mean of $\varepsilon_{is}$.}

\cmtrx{ The adaptivity of the design, with the termination time depending on early values of the outcomes, comes at no cost in the estimation of $\tau$ in the following sense. Suppose we run an adaptive experiment, with a distribution of termination times. Now suppose we compare this to a series of non-adaptive experiments with the same distribution of termination times. This series of non-adaptive experiments is not actually feasible, because it depends on values we do not know {\it ex ante}. Nevertheless, this series of experiments does not do better than our proposed adaptive experiment, in the sense that the average of the variances is the same as that for our proposed adaptive experiment. This result is surprising as the adaptive nature in choosing treatment designs and in experiment termination does not affect the estimation efficiency of $\hat{\tau}_{\mathrm{all},\tilde{T}}$. In Section (ref), we empirically show that the results of Theorem (ref) are valid for a moderate $N$.}

Moreover, as our adaptive treatment decisions seek to maximize $g_\tau(\cdot)$, we expect $g_\tau(\bm{\omega}_{\mathrm{all},1:\tilde{T}},\tilde{T}) > g_\tau(\bm{\omega}_{\mathrm{bm},1:\tilde{T}},\tilde{T})$ so that $\tau$ can be estimated more efficiently than the benchmark design ({\it i.e.}, $\mathrm{Prec}(\hat{\tau}_{\mathrm{all},\tilde{T}}) > ({N\tilde{T}}/{\sigma^2_\varepsilon}) \cdot g_\tau(\bm{\omega}_{\mathrm{bm},1:\tilde{T}},\tilde{T})$). This is shown in Proposition (ref) in Section (ref) for a large $N$, and is empirically demonstrated in Section (ref) for a moderate $N$. }

Extension to Carryover Model with Covariates

PGAE and Theorem (ref) can be easily extended to the specification with $\ell > 0$ and with $\*X_i$. We can continue using the within estimator for $\tau_0, \cdots, \tau_\ell$ by regressing $\dot{Y}_{it}$ on $\dot{z}_{it}, \cdots, \dot{z}_{i,t-\ell}, \dot{\*X}_i$. For the experiment termination rule, we could generalize (ref) to $\mathrm{tr} \big(\mathrm{Prec}(\hat{\bm{\tau}}) \big) \geq c$ or other criteria based on the objectives discussed in Section (ref). From Lemma (ref) in Section (ref), the only unknown parameter in the termination rule is $\sigma_\varepsilon^2$. Furthermore, we can partition units into strata based on $\*X_i$. For each stratum, we can then continue using PGAE to construct an empirical distribution about $\tilde{T}$, make adaptive treatment decisions, and sequentially decide whether to terminate the experiment. We can use a similar proof to show that results in Section (ref) continue to hold with $\hat\tau$ replaced by $\hat\tau_0, \cdots, \hat\tau_\ell$ and $g_{\tau}(\bm{\omega},\tilde{T}) $ replaced by a matrix depending on $\bm{\omega}$ and $\tilde{T}$ (see Lemma (ref) for its definition).

{\color{Black} If the specification has $\*u_i$, there are multiple approaches to proceed. First, as a simple solution, we can ignore $\*u_i$ and run PGAE as the case without $\*u_i$. This approach can be shown to be valid under suitable assumptions.\footnote{For example, the suitable assumptions can be $\*u_i$ is mean zero and i.i.d. in $i$, and $\*v_s$ is i.i.d. in $s$.} We can further improve the precision of $\hat{\tau}_{\mathrm{all},\tilde{T}}$ by re-estimating $\tau$ using GLS post-experiment.

Second, as a more efficient solution, if we have historical data, then we can use it to estimate $\*u_i$, partition units into strata, and run PGAE on each stratum. With abundant historical data, we can precisely estimate $\sigma_\varepsilon^2$, and therefore the minimum duration to achieve a certain precision threshold. For this case, we do not need an adaptive experiment and can instead run the non-adaptive experiment with the estimated minimum duration.}

\setcounter{equation}{0} {\thesection.\arabic{equation}} \oldsection{Empirical Applications}

We run synthetic experiments on multiple real data sets to study our solutions in Sections (ref) and (ref). \footnote{Our code is available at \url{https://github.com/ruoxuanxiong/staggered_rollout_design}.} First, we describe the data sets that we study from multiple domains in Section (ref). Next, in Section (ref), we show that for non-adaptive experiments, our solutions from Section (ref) require less than 50% of the sample size to achieve the same treatment effect estimation error as the benchmark designs. For adaptive experiments, in Section (ref), we show that our adaptive design from PGAE can improve the precision of treatment effect estimation by more than 20%, on top of the improvements obtained by our non-adaptive designs.

Data Descriptions

Our synthetic experiments are run on four different data sets. The first one, MarketScan Research Databases, is used for the empirical results of this section. As robustness checks, the same results are shown on the remaining three data sets in Section (ref).

\paragraph{MarketScan research databases.} These databases contain inpatient and outpatient claim records. Focusing on influenza as the primary diagnosis, there are 21,277 inpatient admissions versus 9,678,572 outpatient visits in the databases. We denote all of these as influenza visits. Our outcome variables are monthly flu visit occurrence rates per Metropolitan Statistical Area (MSA) and per thousand patients, which is defined as the ratio of the number of influenza visits among all enrolled patients times $1,000$ for the given month in a given MSA. Moreover, our analysis focuses on the flu peak seasons that are defined as October to April of the next year. We focus on the period from October 2007 to April 2015 as the databases have few observations outside this period. This leaves us with a panel of $185$ MSAs over $56$ months. See Section (ref) for more details.

\paragraph{Other data sets.} The three additional data sets, as described in Section (ref), are home medical visits in 61 cities over 144 weeks, grocery store transactions for 7,130 households over 97 weeks, and Lending Club loans for 956 geographic areas in the US over 139 months.

Non-Adaptive Experiments

First, in Section (ref), we discuss the setup of synthetic experiments and evaluation criteria, and then present the results in Section (ref).

Setup

\\ Here, we first introduce the benchmark designs as well as different versions of our solution, depending on the specifications of the estimator. Then we explain how the synthetic experiment and treatment effect are generated, and then discuss the evaluation metrics that are used.

{\bf Treatment designs.} We consider the following treatment designs. Illustrations of these designs in Figure (ref) of Section (ref) can facilitate the reading.

enumerate• Benchmark treatment designs: \begin{enumerate} • $Z_{\mathrm{ff}}$ (fifty-fifty): $Z_{\mathrm{ff}}$ has 50% control and 50% treated units at every time period. More precisely, $Z_{\mathrm{ff}}$ is a rounded solution when starting with $\omega_s = 0$ for all $s$. • $Z_{\mathrm{ba}}$ (before-after): $Z_{\mathrm{ba}}$ has all units in the control state before halftime and all units in the treatment state after halftime. More precisely, $Z_{\mathrm{ba}}$ is a rounded solution that starts with $\omega_s = -1$ for $s < (T+1)/2$ and equals to $\omega_s = 1$ for $s \geq (T+1)/2$. • $Z_{\mathrm{ffba}}$ (fifty-fifty with before-after): $Z_{\mathrm{ffba}}$ has all units in the control state before halftime and half of the units in the treatment state after halftime. That is, $Z_{\mathrm{ffba}}$ is a rounded solution that starts with $\omega_s = -1$ for $s < (T+1)/2$ and has $\omega_s = 0$ for $s \geq (T+1)/2$. $Z_\mathrm{ffba}$ combines $Z_\mathrm{ff}$ and $Z_\mathrm{ba}$, and has the simultaneous treatment adoption pattern. \end{enumerate} • Variations of $Z_\textnormal{opt}$ from Section (ref): In order to assess the benefit of various features of our specification (ref), we consider three different designs. Each one is a variant of $Z_\textnormal{opt}$, but with different specifications of the estimator. \begin{enumerate} • $Z_{\textnormal{opt}, \mathrm{linear}}$: This design is optimal under the specification (ref) with $\ell = 0$ and without covariates ($d_x = d_u = 0$). This design can in fact be considered as the “state-of-the-art” benchmark design since it is analogous to the optimal stepped wedge designs of hemming2015stepped and li2018optimal in which the treated fraction increases linearly in time. We sample $\{A_i\}_{i \in [N]}$ from $\mathbb{A}_0$ defined in (ref), and the sampled $\{A_i\}_{i \in [N]}$ uniquely defines $Z_{\textnormal{opt}, \mathrm{linear}}$. • $Z_{\textnormal{opt}}$: This design a nonlinear staggered design and is optimal under the specification (ref) with $\ell > 0$ and without covariates ($d_x = d_u = 0$). We sample $\{A_i\}_{i \in [N]}$ from $\mathbb{A}_{\ell}$ defined in (ref), and the sampled $\{A_i\}_{i \in [N]}$ uniquely defines $Z_{\textnormal{opt}}$. • $Z_{\textnormal{opt},\mathrm{stratified}}$: This design is also a nonlinear staggered design and is optimal under the specification (ref) with $\ell > 0$ and with discrete-valued latent covariates ($d_x = 0$, $d_u > 0$). The value of this design is only demonstrated when historical control data is available, which is a realistic assumption in practice. \cmtfinal{In our empirical applications, we first estimate $\*u_i$ by singular value decomposition (SVD) using historical data (see “evaluation metrics” below for the construction of historical data). Next we partition units into strata based on estimated $\hat{\*u}_i$ and randomly choose a treatment design that satisfies the conditions in (ref) for each stratum, where the number of strata varies from $2$ to $4$. See Sections (ref) and (ref) for more details.} \end{enumerate}

{\bf Synthetic non-adaptive experimental data.} Since we are not aware of any specific experiment that was performed on the data, we assume the data is the control data ({\it i.e.}, original panel data entries are $Y_{is}(-\bm{1}_{\ell+1})$, for all $i$ and $s$). We then create a hypothetical treatment with instantaneous and lagged effects. Given a treatment design $Z$, the observed outcome (in a hypothetical experiment) for unit $i$ at time $s$ would be (recall that $z_{is} \in \{-1,+1\}$) \[ Y_{is} = Y_{is}(-\bm{1}_{\ell+1}) + \sum_{j = 0}^{\ell} (z_{i,s-j} + 1) \cdot \tau_j\,, \] where $z_{is} = -1$ for $s \leq 0$. For the results presented in Section (ref), $\ell$ is chosen at 2. We consider other values of $\ell$ in Section (ref).

{\bf Evaluation metrics.} Instead of running a single simulation on the entire panel of control data, we select $m$ random sub-blocks of dimension \cmtfinal{$N \times (T_{\mathrm{hist}} + T)$, where the first $T_{\mathrm{hist}}$ periods are historical control data and the synthetic experiment is applied to the last $T$ periods of data}. The estimated $\tau_j$ for $j \in \{0\} \cup [\ell]$ on $k$-th block are denoted as $\hat \tau^{(k)}_{j} $. For each design, we use the non-adaptive experimental data generated by this design to estimate $\tau_0, \cdots, \tau_\ell$ using GLS with specification (ref). As a robustness check, we also compare different estimation methods based on different specifications as well. We report the mean and 95% confidence band of total squared estimation error $\sum_{j = 0}^{\ell} \big( \hat \tau^{(k)}_{j} - \tau_j \big)^2$ which is motivated by the objective of A-optimal design, defined in Section (ref). {\color{Black} Note that for GLS, none of the estimation error metrics depend on the actual value of $\tau_0, \cdots, \tau_\ell$, as discussed in Remark (ref). For illustration purposes, we set the lagged effects to decay linearly in lag, $\tau_{1} = 2 \tau_{0}/3$, $\tau_{2} = \tau_{0}/3$, $\tau_j = 0$ for $j > 2$, and the cumulative effect $|\tau_0 + \tau_1 + \tau_2| = 0.1(NT)^{-1}\sum_{i,t} Y_{it}(-\bm{1}_{\ell+1})$. The latter selection means the cumulative effect has a magnitude that is $10\%$ of the average outcome in the panel. We verify that our results are robust to other values of $\bm{\tau}$ and to a much smaller magnitude of the cumulative effect in Figure (ref). } As a robustness check, we also report the squared estimation error of cumulative effect $\big(\sum_{j = 0}^{\ell} (\hat \tau^{(k)}_{j} - \tau_j) \big)^2$, and metrics related to hypothesis testing, that is, the receiver operating characteristic (ROC) curve and the corresponding area under the curve (AUC) in Section (ref).

Results

\paragraph{Staggered treatment designs outperform benchmark designs.} The left subplot in Figure (ref) shows the total estimation error $\sum_{j}(\hat{\tau}_j - \tau_j)^2$ of our nonlinear staggered design $Z_{\textnormal{opt}}$ and benchmark designs $Z_\mathrm{ff}$, $Z_\mathrm{ba}$ and $Z_\mathrm{ffba}$. Both $Z_{\textnormal{opt}}$ and $Z_\mathrm{ffba}$ consistently and significantly outperform $Z_\mathrm{ba}$ and $Z_\mathrm{ff}$. The design $Z_{\mathrm{ffba}}$, as a combination of $Z_\mathrm{ff}$ and $Z_\mathrm{ba}$, performs significantly better, but is still outperformed by $Z_\textnormal{opt}$. Specifically, by using only 50% of the sample size ($N=25$ versus $N=50$), $Z_\textnormal{opt}$ achieves lower estimation error than $Z_\mathrm{ffba}$.

\paragraph{Nonlinear staggered design outperforms linear staggered design.} The right subplot in Figure (ref) compares our nonlinear staggered design $Z_\textnormal{opt}$ with the linear staggered design $Z_{\textnormal{opt},\mathrm{linear}}$. When $\ell > 0$, $Z_{\textnormal{opt},\mathrm{linear}}$ requires 10% more samples than $Z_\textnormal{opt}$ to achieve the same estimation error. Note that the improvement is solely because of $\ell > 0$. In fact, if $\ell = 0$, the treated fraction of $Z_{\textnormal{opt},\mathrm{linear}}$ is optimal. We show this empirically in Figure (ref) of Section (ref) by observing that $Z_{\textnormal{opt},\mathrm{linear}}$ requires about 5% fewer samples than $Z_\textnormal{opt}$ due to the higher variance of the latter.

figure[figure omitted — 591 chars of source]

\paragraph{Stratification further improves upon the staggered treatment design.} The right subplot in Figure (ref) additionally compares our nonlinear staggered designs without stratification $Z_\textnormal{opt}$ and with stratification $Z_{\textnormal{opt},\mathrm{stratified}}$. Using $Z_{\textnormal{opt},\mathrm{stratified}}$ can further reduce 20% samples to achieve the same total estimation error. Overall, this result suggests the existence of latent covariates in the original data. Therefore, when there are latent covariates, we could use the historical data that contains information about latent covariates to design a stratified experiment.

\paragraph{Robustness to additional data sets.} Figure (ref) in Section (ref) shows that the above three findings continue to hold on the other three data sets, as $N$ is varied. Figure (ref) in Section (ref) shows the above three findings continue to hold on all four data sets, as $T$ is varied.

{\color{Black} \paragraph{Robustness to various specifications of the estimator.}

We compare the performance of various treatment designs, as the model specification varies in Figure (ref) in Section (ref). We show that the specification with $\alpha_i$, $\beta_t$, and $\*u_i$ significantly outperforms the specification where either $\alpha_i$, $\beta_t$, or $\*u_i$ is absent. Moreover, we show that $Z_{\textnormal{opt},\mathrm{stratified}}$ performs best under various specifications. Therefore, both the treatment decisions (design) and specification of the estimator play important roles in reducing the estimation error.}

\paragraph{Robustness to other evaluation metrics.} The above three findings continue to hold when the evaluation metric is the squared estimation error of cumulative effect, as shown in Figure (ref) in Section (ref). Figure (ref) in Section (ref) shows the ROC curve of various designs ({\it i.e.}, power vs. significance level), with AUC reported in Table (ref) in Section (ref). Aligned with other metrics, $Z_{\textnormal{opt},\mathrm{stratified}}$ has consistently higher power than all other designs.

Adaptive Experiments

In this section, we run synthetic adaptive experiments and evaluate adaptive designs produced by PGAE. We describe the experimental setup in Section (ref) and then present the results in Section (ref). We show the finite sample properties of Lemma (ref) in Section (ref). We also show the finite sample properties of Theorem (ref) in Section (ref)\cmtfinal{, which implies the validity of the post-experiment inference using estimates produced by PGAE}.

Setup

\\ Suppose the adaptive experiment can run for a maximum of $T_{\max}$ periods in total and $\ell = 0$.\footnote{The results are robust to the hypothetical intervention with carryover effects, and are available upon request. } The adaptive experiment is terminated if the estimated precision is larger than threshold $c$.

{\bf Treatment designs.} Overall, we consider the following three designs

enumerate• Adaptive design: The design produced by PGAE, with dimension $N \times \tilde{T}$\cmtfinal{, where $\tilde{T}$ is the actual termination time observed in the adaptive experiment}. • Benchmark design: The initial design applied to $\mathcal{S}_{\mathrm{ntu}}$, $\mathcal{S}_{\mathrm{atu},1}$, and $\mathcal{S}_{\mathrm{atu},2}$ with dimension $N \times \tilde{T}$, where for all $s\in[\tilde{T}]$, $N^{-1} \sum_i z_{is} = \omega_{\mathrm{bm},s} = (2s - 1 - T_{\max})/{T_{\max}}$ which is optimal when $\tilde{T} = T_{\max}$ (identical to $Z_{\textnormal{opt}, \mathrm{linear}}$ when $\tilde T = T_{\max}$). • Oracle design: The optimal design for a $\tilde{T}$-period experiment for $\mathcal{S}_{\mathrm{ntu}}$, $\mathcal{S}_{\mathrm{atu},1}$, and $\mathcal{S}_{\mathrm{atu},2}$ with dimension $N \times \tilde{T}$ and $N^{-1} \sum_i z_{is} = ({2s - 1- \tilde{T}})/{\tilde{T}}$ (identical to $Z_{\textnormal{opt}, \mathrm{linear}}$ when $\tilde{T}$ is known ex ante).

Note that the dimensionality of the three designs is the same, so we can make a fair comparison of the performance of these three designs.

{\bf Synthetic adaptive experimental data.} Similar to the synthetic non-adaptive experiments, we assume the original data does not contain any specific treatment that we study. Given a treatment design $Z$, the observed outcome for unit $i$ at time $s$ ($s \leq \tilde{T}$) is $Y_{is} = Y_{is}(-1) + \tau_0 (z_{is} + 1)$.

\paragraph{Evaluation metrics.} As before, we randomly select $m$ blocks, each with dimension $N \times T_{\max}$ from the original control data. We report the mean and 95% confidence band of $\big( \hat{\tau}^{(k)}_{0} - \tau_0 \big)^2$, where $\hat{\tau}^{(k)}_{0}$ is the estimated $\tau_0$ on the synthetic experimental data of dimension $N \times \tilde{T}^{(k)}$ based on the $k$-th block of the original data. $\tilde{T}^{(k)}$ is equal to the value of $\tilde{T}$ for that block; that is, $\tilde{T}$ can vary with $k$.

Results

\\ {\color{Black} We show an empirical distribution of the termination time $\tilde{T}$ in Figure (ref) and the estimation error of $\hat{\tau}_0$ of various designs in Figure (ref). Four observations can be made from these figures.

First, PGAE indeed terminates the experiment early when precision exceeds the threshold. As shown Figure in (ref), when $T_{\max} > 7$, the experiment is always terminated quite early ($\tilde{T}<T_{\max}/2$). This early stopping does not compromise the estimation error as Figure (ref) validates that the stopping rule works correctly and the estimation error of the adaptive design always stays below the variance threshold $1/c$. Looking at the results for different values of the threshold $c$, in Figure (ref) and Figure (ref), we see that the termination time tends to increase with the threshold, which is as expected.

Second, the adaptive design from PGAE consistently reduces the estimation error ({\it i.e.}, improves the precision) compared to the benchmark design ({\it i.e.}, non-adaptive design), where the benchmark design is used as initialization in PGAE. This implies that adaptive treatment decisions in PGAE can be useful in lowering the estimation error post-experiment. The reduction is more substantial for a larger $T_{\max}$. This is because when $T_{\max}$ is larger, the benchmark design is further away from the optimal design; hence there is more room for improvement for the adaptive design. The reduction is more than 20% for $T_{\max} \geq 14$.

Third, the adaptive design consistently has a larger estimation error than the oracle design. The difference between adaptive and oracle designs is primarily due to the loss in precision from not knowing $\tilde{T}$ before the experiment starts. As the benchmark design differs from the oracle design and the treatment decisions are irreversible, the “mistakes” made in early time periods persistently continue to impact later periods. If we seek to narrow down the gap between adaptive and oracle designs, we can increase $N$, so that PGAE can learn the experiment termination time faster and make better treatment decisions early in the experiment.

figure[figure omitted — 851 chars of source]

Finally, we note that there is a different trend in how the estimation error varies with $T_{\max}$ for different designs. For the benchmark design, the estimation error generally increases with $T_{\max}$. This is because as $T_{\max}$ increases, the benchmark design deviates more from the oracle design. For the oracle design, the estimation error consistently decreases with $T_{\max}$. This is because $\tilde{T}$ tends to increase with $T_{\max}$, as shown in Figure (ref), and as a result, the precision of $\hat{\tau}_0$ using the oracle design increases with $\tilde{T}$.\footnote{The precision of $\hat{\tau}_0$ using the oracle design equals to $N (4\tilde{T}^2 - 1)/(3 \tilde T \sigma^2_\varepsilon) $, which increases with $\tilde{T}$.} For the adaptive design, since the algorithm stops as the estimated precision reaches the fixed threshold $c$, we expect the estimation error to generally stay flat for various $T_{\max}$. But since the precision estimation is not exact, and is actually a conservative one, there is no specific pattern for fluctuations in the estimation error.

figure[figure omitted — 664 chars of source]

}

\setcounter{equation}{0} {\thesection.\arabic{equation}} \oldsection{Concluding Remarks} In this paper, we study the optimal design of \cmtfinal{staggered rollout experiments}. \cmtfinal{These experiments} are particularly useful for studying \cmtfinal{the impact of} treatments that have causal effects on both current and future outcomes. Our goal is to optimally make treatment decisions for every unit at every time period, in anticipation of most precisely estimating the average instantaneous and lagged effects. This optimization problem can reduce the sample size requirement and directly minimize the opportunity cost of the experiment in practice. We first study the non-adaptive experiments, where the sample size is fixed and treatment decisions are made pre-experiment. We provide a near-optimal solution to the optimization problem. We further study adaptive experiments, where the experiments can be stopped early if needed. We propose the Precision-Guided Adaptive Experimentation (PGAE) algorithm for adaptive experiments. PGAE makes adaptive treatment decisions and allows for valid post-experiment inference. Finally, synthetic experiments on multiple data sets show that our proposed solutions for non-adaptive and adaptive experiments reduce the opportunity cost of the experiments by over 50%, compared to non-adaptive design benchmarks.

\ECSwitch \ECHead{E-Companion}

longtable[longtable omitted — 2,528 chars of source]

\setcounter{tocdepth}{2}

\setcounter{equation}{0} {\thesection.\arabic{equation}} \oldsection{Supplementary Material for Non-Adaptive Experiments}