Fast Selection from Multiple Treatments: A Sequential Method for Principled Digital Experimentation

1. Introduction

Today, many large digital firms continuously experiment to keep up with competitors and adapt to changes in technology and the mercurial nature of consumer behavior and preference. Firms that use digital experimentation seek to establish causal attribution and risk mitigation for novel product or platform changes (Bojinov and Gupta 2022). Further, empirical evidence has demonstrated performance benefits to technologically advanced firms that adopt thorough experimentation approaches to decision making (e.g., Koning et al. 2022).

The benefits of consistent experimentation are substantial, but participants are not an infinite resource. Even on large platforms with many users, the large number of simultaneous experiments and typically small effect sizes—which require considerable sample sizes to be detected reliably—can put a strain on participant resources. Scarcity of participants compels companies to use strategies to get the most out of a limited resource. One approach is to apply many treatments to each participant simultaneously, typically via an orthogonal test plane. Abbasi et al. (2025) demonstrate the unintended consequences of this approach, namely, that the unaccounted-for interaction effects can bias marginal treatment effect estimates, unbeknownst to product managers or scientific partners. One way to sidestep this issue is to assign individuals to only one treatment when tests for more substantial product/platform changes are being made (and minimal bias matters most), but this naturally worsens scarcity.

Alternatively, adaptive procedures can be used to minimize sampling. For these methods, the sample size assigned to a given treatment or policy change is determined as the sampling and parameter estimation unfold, usually with the objective of sampling less from low-performing treatments or to stop sampling as soon as there is sufficient evidence to make a decision. However, adaptivity has been shown to bias treatment effect estimates (Nie et al. 2018), and the bias generally worsens the more adaptive a method is. Estimation bias hinders post hoc studies such as heterogeneous treatment effect analysis, which has been noted as one of the key advantages of online experiments (Karahanna et al. 2018). Reliable follow-up studies then require new samples, worsening scarcity.

Scarcity is in tension with the problem of powering studies sufficiently. Kohavi et al. (2022) reveal that power analysis for hypothesis testing is poorly understood and often ineptly applied in practice in digital experimentation, even by A/B testing consulting firms. This often results in greatly underpowered studies, increasing the uncertainty in decision making and rendering the expended resources effectively useless. Powering analysis is also dependent on knowledge of generally unknown population parameters, such as the variance of the outcome (Cohen 1988), heightening the uncertainty of the test results.

To balance the opposing challenges of scarcity and power for accurate decision making, we propose that firms should use methods that reliably identify the smallest amount of experimentation necessary for accurate decision making while biasing the samples minimally and enabling valuable post hoc analysis. This article proposes such an approach in the setting of comparison and selection among multiple treatments against an incumbent treatment (business as usual) or baseline performance requirement. Our method frames this problem in the null-hypothesis significance testing (NHST) framework, allowing the experimenter to separately specify the error rate of rejecting the standard when there is in fact no better alternative (Type I error) and the error rate of failing to select an existing better alternative (Type II error), unlike most methods of selection and comparison with a standard, which generally fix a single probability of correct selection for the standard and alternatives. Our method then tailors a sequential sampling procedure to accurately and precisely estimate the ideal sample size that meets exactly the desired error rates given the experimenter’s definition of sufficient improvement to implement a change. This aligns statistical power and practical significance as defined by the experimenter (i.e., the minimum effect size that would warrant policy change). Critically, the sequential estimation of the nuisance parameters facilitates a precise estimate of the ideal sample size for the task, meaning our method over- or undersamples very little relative to existing similar methods.

Additionally, whereas this method samples adaptively, it uses sample mean information sparingly. Consequently, the treatment effect estimates are substantially less biased in comparison with many other competitive methods. This better enables post hoc analysis of the treatment effects, such as the examination of heterogeneous treatment effects.

In summary, our method has the following desirable and practical properties, making it uniquely suited for the current state of digital experimentation:

1. The method samples minimally, producing results close to the exact power and size specified, where many methods far oversample and overdeliver on the correct selection rate.

2. The necessary sample size is precisely estimated, resulting in a far narrower upper tail in the sample size distribution in comparison with other methods.

3. The method makes minimal distributional assumptions on the outcome.

4. The method biases treatment effect estimates minimally relative to other methods.

5. No hyperparameter tuning: experimenters only supply parameters following from their problem and preferences; other design parameters are optimized by a procedure we provide.

The organization of this article is as follows. In Section 2, we discuss state-of-the-art methods in digital experimentation, especially for best treatment selection. In Section 3, we formalize the problem and propose the sequential best treatment selection method. Section 4 presents a novel optimization approach that determines the derived design parameters for our method, tailoring the sequential method to the experimenter’s preferences on error rates and practical significance thresholds. Section 5 presents the results of simulation studies demonstrating its performance characteristics and comparisons with existing methods on several criteria, and includes an analysis of data from a real online experiment. In Section 6, we give concluding remarks.

2. Relevant Literature

Digital experimentation methods have become increasingly relevant to the management and information systems (IS) literature. These developments generally approach managerial problems with interesting statistical solutions. Some recent examples include (a) handling interference and heterogeneous effects with switchback experiments (Bojinov et al. 2023), (b) optimizing the design of staggered rollouts for estimating both short-term and long-term treatment effects (Xiong et al. 2024), and (c) always-valid inference methods that protect Type I error rates from “A/B peeking” (data-dependent test stopping) (Johari et al. 2021).

This paper proposes a method in a similar vein to statistical methods for management problems. We develop a novel application of NHST in an algorithmic approach for selecting the best of many treatment arms with comparisons to a standard. In the remainder of this section, we review NHST’s use in digital experimentation and survey existing selection methods, to which we will compare our method.

2.1. Null-Hypothesis Significance Testing in Digital Experiments

In contrast with advances in digital experimentation methods, much of the experimentation in practice is statistically or practically naive, in part because of the sometimes-onerous implementation of theoretically powerful methods. This appears to be especially so for sample size planning for statistical goals (Kohavi et al. 2022). Often, power is not directly considered, and a rule of thumb is used, such as setting data collection to conclude after a fixed period of time, such as two weeks (Kohavi et al. 2013). Alternatively, a fixed n is determined in advance using power analysis, which presupposes some (typically unknown) information about the population. Both approaches regularly result in studies with unnecessarily high power or unacceptably low power with respect to the managerial goals of the evaluation (see Maxwell et al. (2008) for a review of sample size planning methods). Kohavi et al. (2022) emphasize the prevalence of underpowered studies and demonstrate the unreliability of the results produced. Thus, they recommend that experimenters use approaches that are high powered for the typically small effects seen in online settings.

A secondary problem in digital NHST practice identified in the IS and operations literatures is Type I error inflation because of “peeking” or data-dependent test stopping (Johari et al. 2017). The evidence of this kind of error inflation in the field is mixed, with some evidence that it is fairly prevalent (Berman et al. 2018) and some evidence that it is less common in business than in academic research (Miller and Hosanagar 2025). To the degree that this is a real issue, studies that are planned for sufficient power can address it by controlling Type I errors.

Despite its wide use in the field, Feit and Berman (2019) contend that NHST is not ideal for the digital experimentation setting and that less robust testing is more cost-effective. They cite the relatively low consequence of Type I errors, NHST’s high sample size requirement, and the lack of a clear decision when the null is not rejected. Whereas their alternative approach seems quite suitable for marginal changes and tweaks, we believe that many online experimentation settings do justify the NHST framework, especially in cases where more substantial changes are proposed (e.g., new recommender algorithms, large user interface (UI) changes) or there are high-stakes outcomes.

Through qualitative analysis, Koning et al. (2022) discover that firms that successfully implement digital experimentation test both incremental adjustments as well as larger product and website changes. When testing these more substantial changes, Type I errors are not so low cost and can, in fact, result in significant losses or missed opportunities if the implementation of a product has negative effects that are not detected in testing. Additionally, when an experiment is appropriately powered, it is clearer that adhering to the incumbent product or user experience is a reasonable decision when the null is not rejected by the test, so the decision under the null is clear.

2.1.1. Sample Size Planning and Sequential Estimation for NHST.

As Kohavi et al. (2022) emphasize, well-powered studies are key to accurate decision making. Sample size planning is a central part of this process, in which experimenters determine the necessary sample size required to make accurate decisions, conditional on the desired power, Type I error, population effect size, and variance (or a variance-standardized effect size of interest; see Cohen 1988). Nuisance parameters like the population outcome variances must be specified a priori to preplan the sample size of a study to achieve a desired level of power at a given Type I error rate for a given effect size. Population parameters are often difficult to provide because the population parameters are generally unknown. They are often estimated with a small pilot sample, which will tend to render high-variance estimates. Nelson and Goldsman (2001) and Rollin and Chen (2005) exemplify this approach in the best treatment selection literature.

Sequential analysis circumvents the difficulty of presupposing (or estimating with a small pilot sample) unknowable population values by determining sample size a posteriori instead of a priori: samples are drawn successively until a prespecified condition, called the stopping rule, is satisfied. The procedure can then conclude dynamically (in a data-dependent fashion) while controlling desirable statistical properties, and the stopping rule can use all samples up to the stopping point to estimate nuisance parameters such as the variance, which has benefits such as generating precise estimates of the ideal sample size. Notable examples from selection and other statistical domains include Kim (2005), Bandyopadhyay and Biswas (2015), and Kelley et al. (2018).

2.2. Treatment Selection Methods

The problem of selecting the best from several alternatives or arms is addressed by two related but distinct literatures: ranking and selection (R&S) and pure-exploration multiarmed bandit (MAB) algorithms. Although these two domains make different statistical assumptions about the alternatives and seek somewhat different theoretical results, both domains seek to efficiently select the best treatment at a user-specified probability of correct selection. Here, we discuss existing methods from both domains and contrast them with our approach.

A key element shared by many of these methods is the indifference margin. Practically, distinguishing between two arbitrarily close treatments can require an enormous amount of sampling for little gain. Accordingly, the indifference margin is a minimal distance worth detecting, akin to the effect size of interest used to plan sample size in power analysis. If the difference between the mean of the second-best treatment and the mean of the best treatment is less than the indifference margin, then the second-best treatment may be selected at a higher rate than requested because it is deemed “close enough” to the best. In R&S, this margin is generally referred to as the indifference-zone (IZ), whereas MAB algorithms commonly denote it with $ϵ$ (Mannor and Tsitsiklis 2004, Even-Dar et al. 2006).

Without specifying a minimal effect difference worth detecting, a selection algorithm’s goal will be to distinguish between arbitrarily close treatments. Few R&S methods take this approach (see Fan et al. (2016) and related work for a notable exception), whereas many algorithms in the MAB literature do. We believe that this objective is generally not relevant in the digital experimentation setting because the cost of distinguishing between arms with arbitrarily small differences can be extremely high and the returns are as small as the differences.

2.2.1. Ranking and Selection Literature.

Our approach is closely related to a family of R&S algorithms termed “comparisons with a standard,” where a set of alternative proposed treatments is sampled to determine whether any are superior to a standard or existing treatment. For example, Nelson and Goldsman (2001) studied this problem and proposed a two-stage method to estimate the sample size necessary to distinguish the best from a standard and the inferior treatments, akin to fixed sample size planning for a selection task. This two-step approach can be improved with sequential sampling, where the ideal sample size is estimated as samples sequentially arrive, using far more information to precisely determine the total number of samples needed.

Kim (2005) developed a sequential procedure for selecting the best treatment among alternatives and a standard. Her approach uses the initial sample size to estimate the variances of each arm and sequentially samples from each arm, adaptively dropping out apparently inferior treatments. Sequential sampling approaches to R&S can update estimates of nuisance parameters (such as outcome variances) as they go, allowing stopping decisions to be based on high-precision estimates. This will generally reduce the variation in stopping, reducing drastic over- or undersampling because of noisy nuisance parameter estimates. However, this approach requires additional theoretical care. Both the KN++ algorithm proposed by Kim and Nelson (2006) and our method presented in this paper update sample variances sequentially, allowing for more controlled, lower-variance stopping.

We note that the aforementioned R&S algorithms, as well as ours, are oriented in the frequentist perspective to statistical inference. Related work addresses the problem in a Bayesian framework, which can be advantageous, depending on the context and goals (see, for example, Chick and Frazier 2012, Xie and Frazier 2013). Bayesian approaches tend to optimize the sampling and decision process for the expected value across a distribution of possible realities, whereas the frequentist approach tends to optimize for a worst case (typically specified by the IZ). We believe that specifying the minimum worthwhile difference to detect is more intuitive and easier to apply in the context of digital experimentation for managers than specifying a multidimensional prior distribution.

2.2.2. Multiarmed Bandit Literature.

The selection problem this paper addresses is similar to multiarmed bandit problems (Schwartz et al. 2017, Hu et al. 2022). A subclass of MAB problems geared toward pure exploration, often called best arm identification (BAI) methods, comprises methods specifically designed for the selection of the best treatment alternative and is thus most similar to our setting. Especially similar to our work are the “probably approximately correct” (PAC) BAI algorithms that make guarantees of finding an arm within an indifference margin of the true best arm at a specified error rate.

Among others, Even-Dar et al. (2002) and Mannor and Tsitsiklis (2004) pioneered the treatment of this problem in the MAB literature, providing theoretical upper and lower bounds, respectively, on the sample complexity required to identify a near-best arm, with Even-Dar et al. (2006) providing an algorithm that matches the lower bound up to constants. More recent work has further lowered the theoretical sample complexity and developed more efficient algorithms. Two foremost examples are the lil’UCB algorithm (Jamieson et al. 2014) and the Track-and-Stop strategy (Garivier and Kaufmann 2016).

MAB algorithms generally assume Gaussian (or sub-Gaussian or exponential-family) arms with known nuisance parameter values such as arm outcome variance (e.g., Jamieson et al. 2014, Garivier and Kaufmann 2016, Barrier et al. 2022) or else make minimal assumptions except a known bound on the support (e.g., Even-Dar et al. 2006, Karnin et al. 2013), enabling the use of Hoeffding inequalities in the construction of the correct selection guarantees. In contrast, many R&S algorithms make a distributional assumption (e.g., Gaussian) or rely on asymptotics via the central limit theorem (CLT), and many estimate unknown nuisance parameters such as the variance as part of the procedure.

MAB BAI analysis is generally concerned with the lower and upper bounds on the sample complexity, and an optimal algorithm in this setting is one that scales proportionally to the lower bound with respect to the key problem parameters (such as rate of correct selection and indifference margin) “up to constants,” regardless of the size of those constants. This often results in algorithms that are quite conservative in terms of their success rate guarantees, gathering far more samples than their R&S counterparts. To demonstrate this, we include lil’UCB in our main simulation comparisons (Section 5.2) and discuss the sampling behavior of Even-Dar et al. (2006) and Karnin et al. (2013) in our studies with binary outcomes (Section C.3.1 in the Online Appendix and Section 5.4). Although R&S algorithms are, in some dimensions, more practical than MAB approaches, the MAB literature provides strong theoretical analysis of the selection problem, which is invaluable for the field.

3. A Sequential Method for Selecting the Best Treatment

Our procedure builds on a family of two-phase selection procedures, especially Thall et al. (1988), whose selection and testing procedure compares multiple binomial populations, and Rollin and Chen (2005), whose procedure selects the best of several normally distributed populations in the presence of a standard. In these procedures, the first phase screens for the most promising candidate among a set of populations, and the second phase tests whether the chosen population is better than the control population using hypothesis testing. Our method extends the two-phase design developed in these works by implementing a sequential approach and handling the broad class of outcome distributions with finite variances.

Shifting to a sequential approach from the fixed sample size mechanism leverages the sequential arrival of participants in the digital setting and enables conclusion of the experiment with fewer total samples in expectation than nonsequential approaches, as shown in Mukhopadhyay and De Silva (2009, chapter 14) and demonstrated here in simulation studies by comparison with Nelson and Goldsman (2001) and Rollin and Chen (2005) in Section 5 and Online Appendix C.4, respectively. Additionally, sequential sampling enables estimation of the outcome variances with nearly the full final sample, enabling a precise estimate of the final sample size and thus reducing the probability that the procedure over- or undersamples by a substantial amount relative to approaches that rely on pilot samples (e.g., Nelson and Goldsman 2001, Rollin and Chen 2005) and sequential approaches that do not update sample variances (e.g., Kim 2005).

Our method has two design features that make it more suitable for many digital experimentation problems relative to many existing methods. First, the method gives the standard/business as usual a special status with respect to both error rates and minimal significant difference, allowing for a separate indifference margin for selecting an alternative against the standard and a separate error rate for when the standard should be selected (the Type I error rate in our framework) from the error rate for the alternative arm hypotheses (Type II error). These components allow the experimenter to tailor the sequential procedure to one’s conservatism preferences with respect to the standard.

The second feature is that the treatments are not tested for differences at every step. Instead, inferior alternatives are tested and removed only once, partway through the procedure. More adaptive approaches enable stronger early detection and removal of poor treatments, but this greatly impacts the data, biasing treatment effects and other post hoc analyses (Nie et al. 2018), as we demonstrate in simulations (see Section 5.2.3). This stifles follow-up analysis following the experiment. Higher adaptivity also results in higher sample size variation (see Section 5.2.2).

3.1. Formulation

Suppose we want to compare K treatments or policies $P_1,P_2,\dots ,P_K$ (henceforth, “alternatives”), with outcomes having unknown population means ${\mu }_1,{\mu }_2,\dots ,{\mu }_K$ and unknown, finite population variances ${\sigma }_1^2,{\sigma }_2^2,\dots ,{\sigma }_K^2$ against a baseline with known mean ${\mu }_0$. The baseline mean ${\mu }_0$ may represent the mean performance of an existing policy (business as usual) or a baseline performance standard set by a manager. Denote $\mathit{\theta }$ a vector of $\mu$’s specifying a given “configuration” of alternative means, $\mathit{\theta }={[{\mu }_1,{\mu }_2,\dots ,{\mu }_K]}^{⊺}$, and denote the sorted values of the unknown alternative outcome means by ${\mu }_{(1)}\le {\mu }_{(2)}\le \cdots \le {\mu }_{(K)}$. In our two-phase sequential process, our objective is first to determine which alternative has the greatest mean and whether that candidate shows promise to be superior to the baseline (phase I) and, if such a candidate is identified, learn whether the alternative candidate is a significant improvement over the standard (in phase II). The method controls error rates for the joint performance of these tasks so that the decisions are accurate with respect to the overall decision process. This is accomplished by sequentially sampling to dynamically determine $N_i$, the number of samples for alternative $P_i$ for $i\in \{1,\dots ,K\}$ in phase I, and M, the number of samples of the candidate in phase II, given a candidate is selected in phase I.

The experimenter provides two design parameters, ${\psi }_1$ and ${\psi }_2$, $0\le {\psi }_1<{\psi }_2$, that specify the effect sizes of interest in the tasks of comparison and selection. Specifically, ${\psi }_2$ specifies the margin for distinguishing alternative treatments from the baseline (comparison), whereas the difference ${\psi }_2-{\psi }_1$ specifies the indifference margin desired for distinguishing the best alternative from the inferior alternatives (selection). These values jointly specify the closest that the alternative means can be to one another and to ${\mu }_0$ such that the procedure will select the best treatment and determine it to be superior to ${\mu }_0$ at the rate specified by the experimenter.

Alternative $P_i$ constitutes a practically significant improvement over the standard if ${\mu }_i\ge {\mu }_0+{\psi }_2$. Supplying ${\psi }_2$ allows the experimenter to specify the effect size that is practically meaningful to the decision context. This is the minimum effect for which the experimenter would consider replacing the existing standard with an alternative. For example, ${\psi }_2$ may incorporate switching costs or a risk-aversion heuristic. Any alternative with a practically significant improvement or more is considered an “acceptable” selection. The procedure samples to power for the “worst-case” acceptable improvement of ${\mu }_{(K)}={\mu }_0+{\psi }_2$.

Similar to the role of ${\psi }_2$ as the margin between an acceptable alternative and ${\mu }_0$, the margin ${\psi }_2-{\psi }_1$ serves selection among the alternatives by specifying the width of the IZ. An alternative $P_i$ such that ${\mu }_{(K)}-{\mu }_i\ge {\psi }_2-{\psi }_1$ will be reliably distinguished from the best alternative by the selection process. If the second-best alternative is closer than $({\psi }_2-{\psi }_1)$ to the best alternative, the procedure may not sample sufficiently to select the best alternative at the desired success rate. Intuitively, an experimenter will not want to devote enormous resources to distinguishing between similarly successful alternatives. The experimenter should then set ${\psi }_1$ at a level such that $({\psi }_2-{\psi }_1)$ is the smallest consequential difference that the experimenter is willing to expend resources to detect (we provide further guidance for practice in Online Appendix D.1). Table 1 collects the key notation for our formulation, and Figure 1 visualizes the indifference-zone and practical significance regions.

Note that most “comparisons with a standard” algorithms essentially use the same “effect size” for distinguishing between the alternatives as for distinguishing the alternatives from the standard. In our framework, this is equivalent to setting ${\psi }_1=0$. With ${\psi }_1=0$, the same margin is used for distinguishing the best alternative from inferior alternatives and for distinguishing the best from the standard. Allowing ${\psi }_1$ to be set by the experimenter enables differentiation among alternatives with a different margin than differentiation between the best alternative and the standard. This is especially valuable if the experimenter expects several of the alternatives to be competitive relative to the standard and desires a good chance of selecting the true best.

We aim to select the best treatment from the K alternatives if it is at least ${\psi }_2$ greater than the business-as-usual performance or minimum performance standard, ${\mu }_0$. No alternative should be selected if no alternative outcome means are a practical improvement over ${\mu }_0$. For brevity, we denote acceptable means (${\mu }_i\ge {\mu }_0+{\psi }_2$) as ${\mu }_i^{+}$ and unacceptable means as ${\mu }_i^{-}$. Given this notation, we express the goal in the form of a null and alternative hypothesis below.

Here, the null hypothesis ($H_0$) indicates that no alternative improves on the standard ${\mu }_0$, and we refer to the specific parameter configuration where all means are equal as ${\mathit{\theta }}_0$. The alternative hypothesis ($H_a$) states that there exists an alternative that is practically superior to the standard, and no unacceptable alternative ${\mu }_i^{-}$ is so close to the acceptable alternatives that they are too hard to reliably distinguish. The two conditions of $H_a$ jointly specify that the parameters are oriented in such a way that we will achieve power of at least $1-\beta$. Any $\mathit{\theta }$ that satisfies $H_a$ is considered to be in the preference zone (PZ). The “worst-case” value of $\mathit{\theta }\in PZ$ that minimizes the power to $1-\beta$ is the case where both conditions of $H_a$ are tight. This condition is called the least favorable configuration (LFC) and is defined

For any (all else equal) looser configurations, where ${\mu }_{(K)}$ is strictly greater than ${\mu }_0+{\psi }_2$ or ${\mu }_i^{+}-{\mu }_j^{-}>{\psi }_2-{\psi }_1,\hspace{0.17em}j\ne i$, the test power will be greater than the experimenter-specified level because this is an easier-to-detect case.¹ If, on the other hand, either condition is violated while the other is tight, power will be less than the nominal level.

3.2. Tailoring the Algorithm via Derived Design Parameters

To guarantee decision success at the experimenter-specified rates $1-\alpha$ and $1-\beta$, we tailor our algorithm to the experimenter’s design using a set of derived design parameters. Parameters $h_1$ and $d_1$ control sampling in phase 1 of the algorithm by mapping the sequentially updated variance estimates to an estimate of the sample size that provides sufficient precision for the experimenter’s objectives. This sample size estimate operates as a stopping rule. Likewise, $h_2$ and $d_2$ control sampling in phase II. The d and h parameters are derived conditional on the experimenter-selected design parameters ${\psi }_1$, ${\psi }_2$, $\alpha$ and $\beta$, and K, in accordance with Theorem 1.

Theorem 1.

Let ${\psi }_1$ and ${\psi }_2$ satisfy $0\le {\psi }_1<{\psi }_2$ with ${\psi }_1/{\psi }_2\to c$ for some fixed $c\in [0,1)$ as ${\psi }_2\to 0$. Also, let $\alpha \in (0,1)$ and $\beta \in (0,1)$ denote the desired asymptotic Type I and Type II error probabilities, respectively. Let PZ denote the preference zone. Then, for any $\mathit{\theta }\in PZ$ and as ${\psi }_2\to 0$, the proposed test in Algorithm 1 achieves asymptotic power $1-\beta$ and asymptotic size $\alpha$ provided that the design parameters $h_1$, $h_2$, $d_1$, and $d_2$ are chosen to satisfy

$${\displaystyle {\int }_{h_1(1-d_1)}^{\infty }{\mathrm{\Phi }}^{K-1}}\left(z+\frac{h_1{d}_1({\psi }_2-{\psi }_1)}{{\psi }_2}\right)\mathrm{\Phi }\left(\frac{h_2{d}_2}{h_1{d}_1}z-2h_2+2h_2{d}_2\right)\varphi (z)\hspace{0.17em}dz=1-\beta ,$$(1)

$$K{\displaystyle {\int }_{h_1}^{\infty }{\mathrm{\Phi }}^{K-1}}(z)\mathrm{\Phi }\left(\frac{h_2{d}_2}{h_1{d}_1}z-2h_2\right)\varphi (z)\hspace{0.17em}dz=\alpha ,$$(2)

where $\mathrm{\Phi }(·)$ and $\varphi (·)$ denote the standard normal distribution function and density function.

Essentially, Equation (1) controls the rate of Type II errors, whereas Equation (2) controls the rate of Type I errors. Specifically, the left side of Equation (1) expresses the probability, given ${\mathit{\theta }}_{\mathit{\text{LFC}}}$, that the lone acceptable alternative treatment is selected as superior to the other alternatives and practically superior to the baseline as a function of $d_1$, $d_2$, $h_1$, and $h_2$. Likewise, the left side of Equation (2) expresses the probability that any alternative is deemed acceptable relative to the baseline, given $H_0$ (and thus, no alternative is acceptable). We set each of these quantities equal to the experimenter-specified rates and solve for $d_1$, $d_2$, $h_1$, and $h_2$. Using parameter values that solve this system, the algorithm will sample sufficiently to control the error rates. The proof of Theorem 1 can be found in Online Appendix A. We present our optimization approach to solving this system of equations in Section 4.

Theorem 1 relies on the convergence of the sample means to normality, enabling the use of $\varphi$ and $\mathrm{\Phi }$ in Theorem 1. In Section 3.3, after presenting the algorithm, we introduce Lemma 1, a random central limit theorem ensuring Theorem 1’s validity in our sequential setting. Theorem 1 relies on the theoretical asymptotic framework as ${\psi }_2\to 0$, though ${\psi }_2$ is a fixed value provided by the experimenter, akin to letting $N\to \infty$ to invoke the classical central limit theorem in fixed sample settings. In Section 5, we show the efficacy of our method with realistic values of ${\psi }_2$.

3.3. Sequential Algorithm

In this section, we provide a walkthrough and discussion of the two-phase sequential method. Algorithm 1 shows a high-level overview at a glance.

3.3.1. Phase I: Selection.

Given $d_1$, $d_2$, $h_1$, and $h_2$ satisfying Equations (1) and (2), the ideal sample size for alternative $P_i$ for phase I is

which is unknown because it is a function of the unknown parameter ${\sigma }_i^2$. We first take a pilot sample and then sample sequentially to precisely estimate ${\sigma }_i^2$ as $S_i^2$ and ${\nu }_i$ as $N_i$.

3.3.1.1. Pilot Sample Stage.

Draw $n_0$ observations $X_{i1},\dots ,X_{i{n}_0}$ from each of the K treatments $P_1,P_2,\dots ,P_K$. To facilitate a normal approximation of the sample means, we set a minimum sample size of $n_0=30$. This value serves as a practical lower bound and is used in most of our simulation studies. Larger values should be chosen if the data are far from normally distributed (later, our experiment using binary outcomes uses $n_0=100$ samples, which is still small in most digital experimentation contexts). The initial sample sizes chosen for all the treatments are equal ($n_{0i}=n_0$).

Based on this pilot sample of size $n_0$ from the ith treatment, ${\sigma }_i^2$ is estimated using the sample variance $s_{i{n}_0}^2$. This value is used in an initial estimate of the required sample size for $P_i$.

If $n_0<(d_1^2{h}_1^2/{\psi }_2^2)(s_{i{n}_0}^2+1/n_0)$ for all K alternatives, then proceed to the sequential stage. Otherwise, if $n_0\ge (d_1^2{h}_1^2/{\psi }_2^2)(s_{i{n}_0}^2+1/n_0)$ for any $P_i$, $P_i$ should not be sampled any further in phase I. If $n_0\ge (d_1^2{h}_1^2/{\psi }_2^2)(s_{i{n}_0}^2+1/n_0)$ for all $i=1,2,\dots ,k$, stop sampling completely, set the final phase I sample size equal to ${\sum }_{i=1}^K{n}_0$, and proceed to the decision stage of phase I.

3.3.1.2. Sequential Stage.

Here, sequential sampling begins. The method allocates samples according to each treatment’s estimated variance, $s_{i{n}_0}^2$. Intuitively, a larger number of samples are needed for treatments with larger outcome variances to achieve the necessary precision for selection.

At each sequential step, we collect a fixed number $n^{\prime }(\ge 1)$ of additional observations for each treatment that has not yet stopped sampling. This value of $n^{\prime }$ is specified in advance by the experimenter and remains constant throughout the sampling procedure. It may simply depend on the manner in which participants are arriving and typically can be set to one in digital experimentation settings. Larger values risk oversampling by $n^{\prime }-1$.

Using $n_0+n^{\prime }$ observations from each treatment, we find the corresponding sample variance, $s_{i{n}_0+n^{\prime }}^2$. If $n_0+n^{\prime }\ge (d_1^2{h}_1^2/{\psi }_2^2)(s_{i{n}_0+n^{\prime }}^2+1/(n_0+n^{\prime }))$, stop further sampling of $P_i$. If this condition is met for all alternatives, set the final sample size equal to ${\sum }_{i=1}^K(n_0+n^{\prime })$. If this condition is not met for any $P_i$, then continue the sampling process by sampling $n^{\prime }$ more observations from $P_i$.

This process of collecting $n^{\prime }$ observations continues for $P_i$ until $N_i\ge (d_1^2{h}_1^2/{\psi }_2^2)(s_{i{N}_i}^2+1/N_i)$, where $N_i$ is the total number of samples for $P_i$. At this stage, we stop further sampling and report that the final phase I sample size for $P_i$ is $N_i$.

Based on the algorithm just outlined, a sequential stopping rule for each of the K treatments can be defined as follows:

where the term $1/n_i$ is a correction term that ensures that the sampling process does not stop too early. For details about the correction term, refer to Chattopadhyay and Kelley (2016) or Darku et al. (2020). Note that $(s_{i{n}_i}^2+n_i^{-1})$ converges to ${\sigma }_i^2$ as $n_i\to \infty$.

Once the stopping rule is satisfied, using $N_i$ observations drawn from treatment $P_i$, we compute the sample weighted mean for each of the K treatments:

where

These weighted means and the value of $c_{1i}$ ensure that the variances of the weighted means are the same and constant, which is crucial for Lemma 1 shown at the end of this section, which, in turn, enables the asymptotic guarantee made by Theorem 1. The weights are derived and further discussed in Online Appendix B.

Denote the ordered values of the weighted means of the observations from all the K treatments as ${\overline{X}}_{(1)}<{\overline{X}}_{(2)}\dots <{\overline{X}}_{(K)}$.

3.3.1.3. Decision Stage.

If ${\overline{X}}_{(K)}\ge {\mu }_0+{\psi }_2/d_1$, then select $P_S$, the treatment with sample mean ${\overline{X}}_{(K)}$ as the candidate alternative, and proceed to phase II. Alternatively, if ${\overline{X}}_{(K)}<{\mu }_0+{\psi }_2/d_1$, then terminate the procedure entirely and select no alternative. $H_0$ is not rejected.

In summary, phase I samples sufficiently from each of the alternatives to select a sufficiently superior arm (${\psi }_2-{\psi }_1$ greater than others) with fixed probability. Phase I concludes by performing an initial test comparing the empirically best alternative with the standard so that if the empirically best alternative is not promising, phase II will not proceed, and samples can be saved.

3.3.2. Phase II: Comparison with the Standard.

If an alternative with potential to be superior to the standard is found in phase I, phase II proceeds. Phase II allows additional samples to be gathered on only the empirically best alternative and more precisely test whether it improves over the standard. This helps to enable the separate specification of error rates $\alpha$ and $\beta$ and reduces sampling of inferior treatments.

At the start of phase II, we have $N_S$ observations from treatment $P_S$, drawn during phase I. Similar to phase I, the ideal sample size for phase II is

3.3.2.1. Pilot Sample Stage.

Draw $m_0$ additional observations $Y_1,Y_2,\dots ,Y_{m_0}$ from selected treatment $P_S$. Based on this pilot sample, an estimate of ${\sigma }_S^2$ is obtained by computing sample variance $s_{m_0}^2$. If $m_0<(d_2^2{h}_2^2/{\psi }_2^2)(s_{m_0}^2+1/m_0)$, then proceed to the next stage. Otherwise, if $m_0\ge (d_2^2{h}_2^2/{\psi }_2^2)(s_{m_0}^2+1/m_0)$, skip the sequential sampling stage and begin phase II’s decision stage.

3.3.2.2. Sequential Stage.

Obtain $m^{\prime }$ additional observations from $P_S$ and update the sample variance as $s_{m_0+m^{\prime }}^2$. If $m_0+m^{\prime }\ge (d_2^2{h}_2^2/{\psi }_2^2)(s_{m_0+m^{\prime }}^2+1/(m_0+m^{\prime }))$, stop further sampling. Otherwise, continue sampling $m^{\prime }$ observations at a time from treatment $P_S$ until the following stopping rule is satisfied:

Once the stopping rule is satisfied, using the additional M observations from the selected treatment $P_S$, we compute the sample weighted mean ${\overline{Y}}_M$:

where

Both the weighted average and the value of $c_2$ ensure that the variance of the weighted average is constant in the same manner as the weights $a_{ij}$ from phase I.

3.3.2.3. Decision Stage.

If $({\overline{X}}_{\left(KN\right)}+{\overline{Y}}_M)/2\ge {\mu }_0+{\psi }_2/d_2$, then reject $H_0$ and select $P_S$ as the best treatment.

If $({\overline{X}}_{\left(KN\right)}+{\overline{Y}}_M)/2<{\mu }_0+{\psi }_2/d_2$, then do not reject $H_0$, and do not replace baseline treatment.

The terms ${\psi }_2/d_1$ in phase I and ${\psi }_2/d_2$ in phase II are the critical values for detecting a statistically significant difference between the selected alternative and the known standard. Implicitly, they divide this test into two parts, and if the selected alternative exceeds both critical values, then the null hypothesis is rejected. This split is made so that if no alternatives are promising at the end of phase I, the procedure may be concluded early without proceeding to phase II.

Algorithm 1 shows the method at a glance. Note that in the digital experimentation context, “draw $n^{\prime }$ samples from alternative $P_i$” is equivalent to “assign $n^{\prime }$ participants to the $i^{th}$ treatment.”

Algorithm 1

(Sequential Best Treatment Selection)

Determine Parameters: Choose ${\mu }_0$, ${\psi }_1$, ${\psi }_2$, $\alpha$, $\beta$, $n_0$, $m_0$, $n^{\prime }$, and $m^{\prime }$.

Phase I Selection from alternatives

 $N_i\leftarrow n_0$;

 Pilot Stage:

   Draw $N_i$ samples from each alternative; estimate $s_{N_i}^2$.

   Find the Derived Design Parameters

     $h_1,h_2,d_1,d_2\leftarrow$ dp-select$(s_{N_i}^2,{\mu }_0,{\psi }_1,{\psi }_2,K,\alpha ,\beta )$ (see Section 4).

 Sequential Stage:

   for i in 1: K do

     while $N_i<(d_1^2{h}_1^2/{\psi }_2^2)(s_{N_i}^2+1/N_i)$ do

       Draw $n^{\prime }$ samples from $P_i$;

       $N_i\leftarrow N_i+n^{\prime }$

     end

     ${\overline{X}}_{i{N}_i}\leftarrow {\sum }_{j=1}^{N_i}{a}_{ij}{X}_{ij}$ for all i ($a_{ij}$’s are defined in Equation (5))

   end

 Decision Stage:

   Order estimates: ${\overline{X}}_{(1)}\le {\overline{X}}_{(2)}\le \dots \le {\overline{X}}_{(K)}$;

   if ${\overline{X}}_{(K)}<{\mu }_0+{\psi }_2/d_1$ then

     Fail to reject null hypothesis and end algorithm.

   end

   $P_S\leftarrow$ treatment s.t. ${\overline{X}}_{SN}={\overline{X}}_{(NK)}$

Phase II Test best alternative against standard

 Pilot Stage:

   $M\leftarrow m_0$

   Sample M times from treatment $P_S$

 Sequential Stage:

   while $M<(d_2^2{h}_2^2/{\psi }_2^2)(s_m^2+1/m)$ do

     Draw $m^{\prime }$ samples from $P_S$;

     $M\leftarrow M+m^{\prime }$;

   end

 Decision Stage:

   ${\overline{Y}}_M={\sum }_{j=1}^M{b}_j{Y}_j$ ($b_j$’s are defined in Equation (8))

   if $({\overline{X}}_{(SN)}+{\overline{Y}}_M)/2\ge {\mu }_0+{\psi }_2/d_2$ then

     Reject $H_0$, $P_S$ is the best treatment.

   else

     Fail to reject $H_0$, do not replace incumbent treatment.

   end

3.3.3. Lemma 1: Random Central Limit Theorem.

Statistical guarantees in our setting require special care because the final sample sizes $N_i$ and M are random variables determined by the data through sequentially updated sample variance estimates. In such settings, the classical CLT, which assumes a deterministic sample size, is not directly sufficient. Moreover, the underlying outcome distributions may be unknown and need not be normal. Therefore, to ensure the validity of Theorem 1, we establish a random CLT that characterizes the asymptotic normality of our estimators when the sample sizes are random.

Lemma 1.

Provided that the uniform continuity in probability holds, the standardized version of the estimators ${\overline{X}}_{i{N}_i},i=1,2,\dots ,K$ and ${\overline{Y}}_M$ of respective parameters ${\mu }_i$, and ${\mu }_{(K)}$ are such that

$$Z_i=\frac{({\overline{X}}_{i{N}_i}-{\mu }_i)}{\sqrt{s_i^2\sum _{j=1}^{N_i}{a}_{ij}^2}}\underset{d}{\to }N(0,1) \text{for} i=1,2,\dots ,K \text{and} Z_S=\frac{({\overline{Y}}_M-{\mu }_{(K)})}{\sqrt{s_M^2\sum _{j=1}^M{b}_j^2}}\underset{d}{\to }N(0,1) \text{as} {\psi }_2\to 0.$$

The proof of Lemma 1 can be found in Online Appendix A, along with the proof of Theorem 1, where the role of Lemma 1 is made technically explicit.

3.4. Discussion of the Algorithm

The stopping rules in both phases of the algorithm evaluate whether the current sample size from a treatment exceeds the dynamically estimated threshold required for reliable inference. Importantly, the estimated threshold on the right-hand side of the stopping rule depends on the sample variance calculated up to that stage. As a result, treatments with higher estimated variability require larger sample sizes. This approach reflects the principle of Neyman allocation (see Neyman 1934, Zhao 2023), a classical idea from stratified sampling theory that prescribes allocating more samples to groups with higher variability to reduce total estimation error. In our setting, this principle is implemented by adaptively adjusting the required sample size for each treatment according to its estimated variance, ensuring that noisier treatments are sampled more heavily. This variance-sensitive allocation enhances statistical efficiency by aligning sampling effort with estimation uncertainty. For estimating the population variance ${\sigma }_i^2$, we use the conservative estimator $s_{n_i}^2+1/n_i$, where the $1/n_i$ term serves as a regularization buffer. This correction is particularly important during early sampling stages, when variance estimates can be unstable. The idea parallels confidence-based elimination strategies in the multiarmed bandit literature (e.g., Even-Dar et al. 2006), where uncertainty margins are used to avoid premature decisions. In our context, the $1/n_i$ term inflates the variance estimate slightly to guard against underestimation, helping ensure that no treatment is dropped too early before its performance is assessed with sufficient confidence.

4. Optimizing the Derived Design Parameters: The dp-Select Procedure

In this section, we present a novel optimization approach to map the experimenter’s test preferences to the derived design parameters using the system of equations from Theorem 1. Recall from Section 3.2 that the derived design parameters $d_1$, $d_2$, $h_1$, and $h_2$ are selected such that the Type I error ($\alpha$) and power ($1-\beta$) are met according to the experimenter-selected design parameters ${\psi }_1$, ${\psi }_2$, and K. Section 3.3 showed that ${\mu }_0+{\psi }_2/d_1$ and ${\mu }_0+{\psi }_2/d_2$ are the critical values for phase I and phase II, respectively. The two tests against these critical values jointly determine whether there is a treatment that is statistically superior to the standard. Given these critical values, the terms $h_1$ and $h_2$ adjust the target precision level of each ${\overline{X}}_i$ so that the required error rates are met.

Any combination of parameters $d_1$, $d_2$, $h_1$, and $h_2$ solving the system of equations in Theorem 1 will provide a valid sequential procedure with respect to the error rates. However, solutions to the equations alone do not provide optimality with respect to minimal sample size. Because there are two equations with four unknown values, the problem is underdetermined, and there are infinite solutions that may give different final sample sizes, conditional on the true, unknown variances via their estimates from the data. In order to find the solutions that both provide the desired error rates and minimize the sample size, we minimize an expression of the expected sample size subject to the error-rate system of equations as constraints. This way, for any experiment (parameterized by K, ${\psi }_1$, ${\psi }_2$, $\alpha$, and $\beta$), our algorithm can provide the design parameters $d_1$, $d_2$, $h_1$, and $h_2$, which are optimal with respect to the error rates as well as sample size.

The indeterminacy of the system of equations is the result of evaluating a multiple-decision process based only on the final decision error rates. In choosing $h_1$ and $h_2$ conditional on fixed critical value parameters $d_1$ and $d_2$, the system is effectively allocating precision (or inversely and relatedly, allocating error probability) between the two decision steps, partially determining the amount of information that will be gathered in each phase. Without minimizing the expected sample size with respect to these parameters, a valid solution (with respect to achieving error rates) could “allocate” $99\%$ of the acceptable error probability to the first stage, requiring that the second stage gather a colossal number of samples to nearly ensure the correct decision is made in that phase. This is clearly undesirable, as it results in a greater number of total samples than if the error probability is more evenly distributed between the steps. Here, we address this issue.

The expression for the expected total sample size, $E(N_{\mathit{\text{total}}})$, comes from the stopping rules

where the identity of the selected treatment is unknown. We use the following function, proportional to the sample estimate of the quantity above, and minimize it subject to the error constraints (Equations (1) and (2)) to find optimal design parameters $d_1$, $d_2$, $h_1$, and $h_2$:where ${\widehat{p}}_{n_0,i}$ is the estimated probability that alternative $P_i$ is selected for the second stage of sampling. $s_{n_0,i}^2$ and ${\widehat{p}}_{n_0,i}$ are calculated for each alternative from the initial $n_0$ samples taken at the beginning of the procedure. These initial samples can be drawn prior to determining the derived design parameters. We minimize f subject to the equations in Theorem 1 in Section 3.2 as constraints. We call this operation “design-parameter selection,” or dp-select. The dp-select procedure is included in the open-source R package for our sequential method.

5. Performance Evaluation

We evaluate the performance of our method through a Monte Carlo simulation study. Following a table showing our method’s performance, we compare our method to other methods from the literature under both normally and nonnormally distributed conditions, where our method is competitive or superior in each dimension we evaluate. Lastly, we show our method’s efficacy in solving a problem with real online experiment data. For simulations with at least one acceptable treatment, we evaluate our method and others on P(CS), the probability of correctly selecting the best alternative treatment. Where applicable, each method is calibrated to the true difference between the best alternative and the inferior ones (which would be unknown in practice). This allows us to empirically compare the selection success rate guarantees and sample sizes for each method we compare, which are all relative to some effect-size/indifference margin.

5.1. Simulation Study on the Size and Power of Our Method

In each condition simulated, data were generated either under the null hypothesis (${\mathit{\theta }}_0:{\mu }_1={\mu }_2={\mu }_3=\cdots ={\mu }_0$) or the alternative hypothesis in the least favorable configuration (${\mathit{\theta }}_{\mathit{\text{LFC}}}:{\mu }_0\le {\mu }_{[1]}={\mu }_{[2]}=\cdots ={\mu }_0+{\psi }_1<{\mu }_0+{\psi }_2={\mu }_{[K]}$). Outcomes are normally distributed, and ${\mu }_0=0$ in all cases. $\alpha$ and $1-\beta$ are the nominal size and power of the test, respectively. For each condition, the empirical Type I error rate is the proportion of simulations rejecting $H_0$ when data are generated under ${\mathit{\theta }}_0$, whereas P(CS) is the rate at which $H_0$ is rejected in favor of ${\mu }_{(K)}$ given ${\mathit{\theta }}_{\mathit{\text{LFC}}}$.

Table 2 shows the results of this simulation study. The sample size columns show the average sample sizes across simulations for phase I samples (one for each treatment, $n_1$, $n_2$, and $n_3$) and the average phase II sample size m. The results show that the empirical Type I error rates for each configuration tested fall well within Bradley (1978)’s liberal criterion range of 0.025 to 0.075, and observed P(CS) is similarly near the desired level.

5.2. Performance Comparisons with Normally Distributed Outcomes

To benchmark our approach against existing R&S and bandit algorithms, we first evaluate performance under normally distributed outcomes. Importantly, our procedure does not require normality and only assumes finite variances. In this and the next subsection, we will show that our method is competitive with or superior to existing methods with respect to several important metrics. In other words, we find that any method matching or slightly exceeding our method’s performance with respect to one metric is strongly outperformed by our method on a different metric.

There are two practical features of our method that cannot be compared in simulation. First, nearly all existing methods in either the R&S or MAB domain use only one error rate, $1-P(CS)$, whereas we allow experimenters to specify both Type I and Type II error rates in our framework. For these comparisons, we constrain our method to the limitations of the others by fixing our Type I and Type II error rates to the same value as the $1-P(CS)$ used by the other methods.

Second, of the methods we compare with, even those that specify a “standard,” do not allow experimenters to specify different margins (or minimal effect sizes of interest) for the task of distinguishing among alternatives versus testing against the standard, which we permit by allowing experimenters to specify nonzero values for ${\psi }_1$. We thus constrain ${\psi }_1=0$ for these comparisons.

We evaluate our method next to popular “comparison with a standard” procedures from the R&S literature, specifically Nelson and Goldsman (2001) and Kim (2005), as well as KN++ (Kim and Nelson 2006),² a fully sequential R&S procedure, and lil’UCB, a state-of-the-art selection procedure from the MAB literature (Jamieson et al. 2014). We compare with lil’UCB because it makes a sub-Gaussian distributional assumption and so is suitable in this setting. Many popular MAB approaches (e.g., Even-Dar et al. 2006, Karnin et al. 2013) assume that rewards are bounded in a known interval in order to apply Hoeffding-type concentration inequalities. Such assumptions are restrictive in settings where outcomes are modeled using unbounded distributions such as the normal distribution. In Section 5.4, we discuss how these MAB methods perform relative to our method in the binary outcome case, where their models are better suited.

Each simulation is run with $K=3$ alternatives. We vary two key parameters in these simulations: the efficacy of the best treatment and experimenter-specified P(CS). As the efficacy of the best treatment increases, it is easier to distinguish from the other treatments, and so the algorithms should sample less. As P(CS) increases, the required precision increases and the sample size with it. We hold the inferior means and standard ${\mu }_0$ fixed at zero, and the indifference margins (where applicable) are set equal to the best treatment’s lift from zero.

5.2.1. Sample Size.

First, we compare the methods in terms of the number of samples used to make a decision. We additionally show the empirical P(CS) for each to evaluate the success rates.

Importantly, our method, like many in the R&S literature, seeks to succeed with a rate near the desired rates for a given orientation of the true means (in our case, specified by the null hypothesis and the LFC for $1-\alpha$ and $1-\beta$, respectively). For these methods, good performance means having an error rate near the desired error rates. In contrast, some methods sample far more in order to nearly ensure that the error rates will be below the desired error rate. From the NHST point of view, this is wasteful and unnecessary sampling, in a sense “buying” more than the power for which the experimenter intends to budget.

Figure 2 transforms the y-axis by ${\text{log}}_{10}$ in order to visually distinguish between the R&S algorithms while including lil’UCB in the chart. The y-axis of Figure 3 shows the difference between the empirical success rate and the P(CS) specified by the experimenter.

Figure 2 shows that our method is competitive in terms of sample size, and Figure 3 shows that our success rate meets the specified P(CS) as ${\psi }_2\to 0$, in accordance with Theorem 1. As ${\psi }_2$ grows, our observed success rate dips slightly, partly because of bias from sequential updates to the variance, where earlier-than-optimal stopping can occur because of chance errors in the sample variance estimate. However, the success rate stays within Bradley’s liberal criterion of a 2.5% range around the target rate (Bradley 1978). This shows that our approach deviates little from the error rate guarantees, even in the nonasymptotic condition.

Among the R&S algorithms, Kim (2005) tends to sample slightly less than ours, and KN++ (Kim and Nelson 2006) samples around the same amount, whereas Nelson and Goldsman (2001) consistently sample somewhat more. Despite sampling slightly less than us, we argue that our superior performance in other dimensions can make our approach more practical and desirable than Kim’s for many digital experimentation applications. We discuss these metrics next.

5.2.2. Sample Size Variation.

Sample size standard deviation is a critical metric, as it measures an algorithm’s tendency to over- or undersample with respect to the objective. Many methods perform well on average but suffer from large upper-tail behavior, regularly sampling far more than necessary. We plot the standard deviation of the sample size (Figure 4) and the $95\%$ quantile across simulations to represent the upper tail (Figure 5). Note that the y-axis here is also transformed using the ${\text{log}}_{10}$ scale. These graphs show that our approach achieves far smaller variation in sample size relative to others.

5.2.3. Effect Estimation Bias.

As Nie et al. (2018) demonstrate, adaptive sampling procedures bias the sampled data. Bias in the samples reduces the reliability of the data for post hoc analyses, like heterogeneous effect discovery, worsening data scarcity. However, not all sampling procedures are biased in the same way, and some bias the data far less than others. Here, we compare the methods in terms of the bias in the average outcome estimates for the best treatment and an inferior alternative treatment. Figure 6 and Figure 7 show the estimation bias results.

Usually, bias is introduced when an algorithm uses sample mean information to make sampling decisions. Nelson and Goldsman (2001) never use sample means to decide which arm to sample or how much. Thus, their sample mean estimates are unbiased, and their performance serves as a benchmark. Kim (2005), KN++ (Kim and Nelson 2006), and lil’UCB (Jamieson et al. 2014) are highly adaptive procedures, using the sample mean information at every step to decide whether to continue sampling. This is beneficial when the primary objective is to stop early, but can greatly bias the sample. Our approach takes a middle ground, where the sample mean information is used once at the conclusion of each phase. This reduces how much the inferior treatments are sampled relative to Nelson and Goldsman (2001), while taking on less bias than highly adaptive algorithms such as Kim (2005) and Jamieson et al. (2014).

Kim (2005)’s method exhibits substantial bias for all arms—substantial positive bias for the acceptable arm and substantial negative bias for the unacceptable ones (compounding the bias on the difference in treatment outcomes, which is often of interest in post hoc analysis). KN++ shows bias in some cases but is generally reduced relative to Kim (2005). In comparison, our method shows a relatively small amount of negative bias for all arms.

The difference in the direction of bias of the best arm between our method and the others is also because of adaptivity differences. Kim, KN++, and lil’UCB will stop sampling at any time when one treatment mean is substantially greater than the others. This results in many cases where sampling ends when the best treatment’s sample mean was biased positively because cases where the sample mean was not biased positively lead to more samples, reducing bias. In contrast, our bias largely comes from the cases where the best treatment mean is not selected as the best at the end of phase I, occurring because of negative errors in estimation in phase I. This results in a small portion of samples with negative bias, which, in turn, bias the sample means across simulation runs.

5.3. Performance Comparisons with Non–Normally Distributed Outcomes

Our method requires only mild moment conditions and does not assume a specific outcome distribution. As a result, it can be applied in a wide range of experimental settings common in digital experimentation, including binary outcomes. Some of the methods we compare against require more stringent assumptions, and their performance is expected to suffer in very nonnormal conditions (e.g., Nelson and Goldsman (2001) and Kim (2005) assume normality). In contrast, methods making a more general distributional assumption should be unaffected, such as Jamieson et al. (2014), which makes a sub-Gaussian assumption on the outcomes, and KN++ (Kim and Nelson 2006), which uses asymptotics and an adapted CLT, akin to our approach. Figure 8 shows P(CS) results for Bernoulli outcomes. We provide further empirical comparisons with these methods in conditions with gamma- and Bernoulli-distributed outcomes in Online Appendix C.3. In general, all methods perform reasonably well in the gamma-distributed condition, whereas only our method and KN++ are able to maintain their statistical guarantees in the Bernoulli-distributed conditions.

5.4. Online Marketing Experiment Demonstration

We partnered with a marketing team that was testing five different ways of advertising an advanced business school degree program. Each advertisement was displayed to several thousand individuals on the web, and the number of individuals who clicked an ad was recorded. Each event where an individual sees an ad is an “impression.” The marketing team sought to determine which ad generated the highest click rate $\left(\#\text{clicks/}\#\text{impressions}\right)$. All click rates (and differences between them) are small, as is typical in digital experiments. Because the marketing campaign was completed before our partnership, we treat the observed data and its click rates as populations from which to sample. We sampled from the click outcomes for each treatment to simulate a sequential experiment on binary outcomes and select the best of the five ads. We repeated the experiment 10,000 times per condition to estimate the mean and variance of the final sample size and observed P(CS).

As discussed previously, for Bernoulli distributions with small p as in this data set, the sample mean converges in distribution to the normal distribution at a slower rate, so we use $n_0,m_0=100$ samples. Note that this is far below the number of samples taken per arm in the first round of MAB methods, such as the median elimination algorithm of Even-Dar et al. (2006) for the small effect sizes in this context.

We suppose that the marketing team has a minimum performance expectation (derived from the break-even expectations balancing ad cost and click value) of one click in 200 impressions, or a click rate of 0.005. They then use this as a baseline performance expectation, which treatment alternatives must dominate (${\mu }_0=0.005$). We then power to detect a superior (acceptable) ad achieving 1.2 clicks out of 100 impressions (${\psi }_2=0.012-0.005=0.007$), they are willing to budget to distinguish between alternatives as close as 0.004 clicks/impressions (${\psi }_1=0.003$),³ and they seek a correct selection rate of 0.9 ($\alpha =\beta =0.1$).

Table 1. Summary of Design Parameters

Table 2. Error Rate and Sample Size Estimates for Our Method in Simulation

Table 3 shows the results of these experiments. For this demonstration, we assume the empirically best treatment in the full data set is the true best, so “success” is defined as selecting that treatment. In the case of binary rather than normal outcomes, we observe a higher rate of correct selection than desired, which means that our sample size is somewhat larger than is necessary, although this is because of, at least in part, the “population” p values being in a position that is slightly easier to detect than the LFC (because the inferior ps are not all exactly equal to ${\mu }_0+{\psi }_1$; some are smaller). Contrast this with the median elimination (Even-Dar et al. 2006) and exponential-gap (Karnin et al. 2013) algorithms: for an experiment with desired $\text{P}(\text{CS})=0.9$ and indifference margin of 0.004, median elimination samples $1/{(0.004/8)}^2·\text{log}(6/0.1)=16,377,378$ (using the natural log) observations from each treatment in the first round alone. The exponential-gap algorithm uses median elimination as a subroutine, and by the third iteration, exponential gap generally eliminates no alternatives (including the standard), and the first round of the median elimination subroutine draws nearly three million samples per arm.

Table 3. Results from Analysis of Digital Ad Experiment Data

Notably, the original experiment consisted of 123,840 impressions across all treatments. Thus, with these small effect sizes and stringent error rate requirements, we sample substantially less than the original experiment.

6. Conclusion

Our sequential method for determining the best treatment is a pragmatic approach that maps closely to a managerial need in digital experimentation: well-powered experiments when participants are scarce. This approach is designed to use near the ideal sample size, the exact number of samples implied by an experimenter’s objectives with respect to error rates ($\alpha$ and $\beta$) and effect margins (${\psi }_2$ and ${\psi }_2-{\psi }_1$), and the sequential estimation of the necessary sample size results in substantially less over- or undersampling than competitor methods. Unlike many other methods in this space, our method requires no hyperparameter tuning: the design parameters correspond directly to the experimenter’s preferences and objectives, whereas the derived design parameters ($h_1$, $h_2$, $d_1$, and $d_2$) are provided by our optimization routine dp-select. Our sequential selection method and procedure for finding the design parameters are implemented in an R package.⁴

The Monte Carlo simulation studies demonstrate that our method is competitive with existing approaches in terms of sample size and superior in terms of sample size precision. The simulations also show that estimation bias is less than Kim (2005) and Kim and Nelson (2006), the comparison methods that are most competitive in terms of sample size. Monte Carlo studies also validate that our method is just as effective in nonnormally distributed conditions when starting sample sizes ($n_0$, $m_0$) are large enough to engage the central limit theorem, whereas some of our competitors do not adapt as easily. Simulations with gamma-distributed populations show resilience to skewness, and simulations with small-p binary data and a study based on data from a real digital experiment show effectiveness given small-p binary outcomes.

Despite its applicability, stopping rule–based methods such as ours can produce practical challenges for experimenters. A primary challenge is that it is possible that the stopping rule is prohibitively out of reach—far beyond the sample budget of the experimenter. This is essentially the case for all R&S and MAB methods we discussed. However, the simplicity of our stopping rule allows experimenters to estimate the remaining sample size needed, and this feature is implemented in our R package. Once the sample size is high enough to have a reasonably precise estimate of the sample variance, it is unlikely that the estimated sample size will change dramatically. It is up to the experimenter’s discretion to stop early in this case. If stopping is done solely based on the projected sample size being well estimated at a prohibitively large value, the samples are not contaminated. In other words, early stopping in this case does not prohibit the further usage of those samples, although the result of our method must be a failure to reject the null hypothesis in that case.

One limitation of our approach relative to the highly adaptive sampling methods such as Kim (2005), Kim and Nelson (2006), and Jamieson et al. (2014) is that we will generally not stop poorly performing treatments as early. This is a direct trade-off with the low sample size variance and low effect estimate bias, both of which generally increase with adaptivity. If an experimenter has a prior belief that most of the treatments are poor, one of these may be a preferable approach. However, in digital experimentation, very small differences between treatments are quite common, as evidenced by the digital experiment field data analyzed in Section 5.4.

Additionally, our method assumes a known baseline treatment, representing a known business-as-usual mean outcome or baseline performance expectation. In some applications, an experimenter may have no baseline performance expectation and have only a noisy estimate of the business-as-usual condition. Adapting our method to use an estimated baseline would make a valuable extension of our work that would require a generalized form of Theorem 1 for such conditions. In lieu of this advancement, experimenters might use an optimistic estimate or upper confidence interval value as ${\mu }_0$ to form a conservative test.

Although we propose a particular method in this work, the managerial problem of balancing scarcity and test power is broader than the problem of selection among alternatives. Experimenters employ a broad set of tests, from simple A/B tests to customer subgroup heterogeneity analyses. Future method development in these domains should increasingly focus on precisely balancing power and resource usage.