# Anytime-Valid Confidence Sequence-based Experimentation Service for Online
A/B Tests

Akash V. Maharaj1, Ritwik Sinha2, David Arbour2, Moumita Sinha1, Simon Liu1 ,
Raghavendra Addanki2, Ian Waudby-Smith3, Aaditya Ramdas3 and  Manas Garg1,
Viswanathan Swaminathan2 1 Adobe Inc.2 Adobe Research3 Carnegie Mellon
University

(2023)

###### Abstract.

A/B testing is the gold standard for testing features to optimize digital
experiences. The ubiquity and ease of conducting A/B tests requires designing
and implementing solutions to safely conduct A/B tests that are readily used
by practitioners (without the need for oversight by trained statisticians). It
is widely known that allowing the experimenter the ability to ”peek” (evaluate
whether the test has reached conclusion) in traditional A/B testing mechanisms
can lead to inflated type-I error and erroneous conclusions. We have released
an experimentation service on the Adobe Experience Platform based on Anytime
Valid Confidence Sequences, that allows continuous monitoring of the A/B test
which can be concluded as soon as the stopping criterion are met. Our
Experimentation Service is full featured, allowing for sample size
calculations, analyzing multiple treatments (A/B/n tests), and the ability to
analyze alternate test statistics like the “lift”. On simulated data, we show
the desirable properties of using the Anytime Valid Methods compared to
alternate sequential methods as well as fixed horizon test statistics. We
further show the properties of our test when compared with fixed horizon
methods on thousands on real experiments conducted in our product.

A/B testing, randomized experiments, p-values, anytime valid, always valid

††copyright: acmcopyright††journalyear: 2023††doi:
XXXXXXX.XXXXXXX††conference: ; Apr 30–May 4, 2023; Austin, TX††price:
15.00††isbn: 978-1-4503-XXXX-X/18/06††ccs: Computer systems organization
Embedded systems††ccs: Computer systems organization Redundancy††ccs: Computer
systems organization Robotics††ccs: Networks Network reliability Figure 1.
Image of an experiment being analyzed with Anytime-Valid Confidence Sequences
in an online A/B test conducted using the Experimentation Service on Adobe
Experience Platform.

##  1\. Introduction

A/B testing — also known as randomized experimentation — is a ubiquitous tool
for testing features to optimize customer experiences on websites, mobile
applications, and in emails. Over the years, there has been a steady and
substantial growth in the use of A/B tests for improving every aspect of a
customer’s experience. This growth has led to systems and solutions that
enable easy implementation of A/B tests on a wide variety of platforms. This
ease of use means that today, developers, product managers, marketers,
analysts, designers and business leaders have access to tools with which then
can conduct A/B tests. While this ubiquity has led to the growth of data-
driven decision-making within organizations, it also means that practitioners
with limited understanding of the nuances of statistical inference have access
to powerful tools that could lead them astray if used improperly.

This democratization of A/B tests motivates the need to provide testing
solutions that can protect practitioners from drawing flawed conclusions from
tests. One aspect of A/B tests that can often lead to flawed conclusions is
“peeking”, where the person conducting the test has multiple opportunities to
compare $p$-values (or confidence intervals) to some threshold (e.g. 0.05).
This is well known to severely inflate the type-I error for so-called “fixed
horizon” A/B tests (johari2015can; deng2016continuous; johari2017peeking;
johari2022always; howard2021time). A fixed horizon A/B test (the traditional
approach to conducting a hypothesis test) is one where the study is designed
by pre-specifying the null and alternate hypotheses, and a desired sample size
is calculated for the test (based on a minimum detectable effect and desired
statistical power, which depend on unknown quantities). Only _when this pre-
specified sample size is reached_ , can the experimenter compare the $p$-value
(or confidence interval) to the appropriate threshold. Furthermore, this
comparison should _only be done once_. Performing these comparisons multiple
times before reaching this sample size, or continuing to collect more data and
performing additional comparisons will drastically inflate the type-I error
(wald1945sequential; robbins1970statistical; johari2015can;
deng2016continuous; johari2017peeking; howard2021time). For randomized
clinical trials — another common application area for randomized experiments —
the Food and Drug Administration (FDA) strongly recommends prespecification of
designs and controlling the number of peeks(fdaguidance), with use of group-
sequential methods (jennison1999group) prominently highlighted. In digital
experimentation, however, most tests do not have a regulatory body enforcing
design pre-specification, and even if they did, the flexibility provided by
the methods we present here — so-called “confidence sequences” — may be
preferable.

The issue of peeking has received attention in web experimentation; in
particular, johari2017peeking studied its adverse consequences, and proposed
solutions based on sequential hypothesis testing and confidence sequences (CS)
(darling1967confidence; robbins1970statistical; howard2021time). However, most
of the literature on CSs focuses on nonasymptotic methods, which have three
major disadvantages even for fixed-horizon settings: (a) they require strong
assumptions, such as a parametric model or other substantial information about
the distribution (robbins1970statistical; howard2021time;
waudby2020estimating; waudby2020confidence), (b) they are typically wider than
asymptotic ones based on $t$-tests or the central limit theorem, and (c) they
take different forms for different problems, whereas asymptotic methods
typically have a universal and closed-form (trivial-to-compute) expression.
Indeed, for the fixed-horizon regime, CLT-based confidence intervals avoid
issues (a)–(c); a natural question is: is there an anytime-valid analogue of
these methods? waudby2021doubly answer this question in the affirmative by
developing “asymptotic confidence sequences” (AsympCS) which we leverage to
build and deliver a simple, general, and efficient experimentation ecosystem
on the Adobe Experience Platform.

Our contributions include the following: 1) we adapt Anytime Valid Confidence
Sequence for the average treatment effect of an A/B test to a scenario where
propensities are not recorded during decision making (leveraging
(waudby2021doubly)) [I am confused here — waudby2021doubly already does this,
and this has been the case since v1 of the paper. What is different here?
–Ian][key difference is that our analytics systems don’t capture propensities,
and can’t (easily) compute the variance of the IPW estimator. So we just
trivially outlined how to calculate the variance using only means and
variances of individual arms –Akash] 2) we propose correcting for multiple
treatment comparisons using the Benjamini-Hochberg procedure
(benjamini1995controlling), [I think you should be using the e-BY paper
instead, to correct for selecting promising CIs, controlling both FCR and FDR:
https://arxiv.org/abs/2203.12572 — it applies to CSs, even with dependence,
etc. –Aaditya] [Just to add to Aaditya’s comment: I believe that moving from
BY to e-BY is not complicated at all (the procedure is almost exactly the
same) so this really is worth it. In particular, the FCR guarantee is valid
under _arbitrary_ dependence, and hence you can even use the same users as
subjects across various A/B tests and still have valid FCR control — this is
_not_ possible with classical BH/BY without additional assumptions /
modifications –Ian][Acknowledged. At time we were implementing we weren’t
aware of this work, hence did not do it. Given pace of things, may take us a
few months to implement this –Akash] 3) we propose a way of computing AsympCS
for alternate test statistics like the “lift” (relative improvement), 4) we
propose a sample size calculator for AsympCS based tests, 5) using
simulations, we show that AsympCS has attractive properties when compared to
alternate strategies commonly used in the industry, 6) on thousands of real
A/B tests conducted within the product, we contrast the advantages of using
AsympCS over fixed horizon tests.

##  2\. Related Work

Fixed horizon tests (Chapter 10 of (wasserman2004all)) involve 1) specifying
the null and alternate hypothesis, 2) the desired type-I error ($\alpha$,
often 5%), 3) a minimum detectable effect, 4) these lead to the sample size
needed to achieve a certain statistical power ($1-\beta$), 5) after waiting to
observe this sample size, the experimenter can compare the test statistic
using the $p$-value (or confidence interval) exactly once. Allowing multiple
peeks leads to increased type-I error ($>\alpha$) (johari2015can). This
challenge has been long known (pocock1977group) and the traditional approach
to addressing it is to strictly define and monitor the times at which the
comparisons are made (pocock1977group; o1979multiple), known as Group
Sequential Methods. Lan and De Mets (LDM) (gordon1983discrete) proposed an
extension of a Group Sequential method that can provides more flexibility. LDM
does require a sample size to be specified, but the times of analysis can be
flexibly determined. Decision boundaries are computed by considering the times
of the prior looks, and the amount of type-I error spent so far.

Recent work has focused on methods that do not require pre-specifying a sample
size or keeping track of the times of these looks. A particularly impactful
contribution in the technology industry was the proposal of a mixture
sequential probability ratio test (mSPRT)(johari2017peeking). This method is
nonasymptotically valid for parametric exponential families, and shows
excellent properties in practice. However its reliance on this parametric
model means that its validity depends on the correctness of the underlying
model. Subsequent work has focused on using Martingale results to define
Anytime Valid Confidence Sequences [Johari’s results also use martingale
techniques and they explicitly derive confidence sequences. Moreover, the
mSPRT-type confidence sequence already appeared in lots of early (60s/70s)
work by Robbins et al. –Ian] (howard2021time; waudby2020confidence;
waudby2020estimating). Formally, a Confidence Sequence is a sequence
$\bar{C}_{n}$ of confidence bounds, such that, $\mathbb{P}(\forall
n\in\mathbb{N}^{+},\theta\in\bar{C}_{n})>(1-\alpha)$. Initial approaches
primarily introduced non-asymptotic confidence sequences. However, these
require strong assumptions or produce wide intervals. The most recent
developments have introduced _asymptotic_ confidence sequences
(waudby2021doubly). AsympCSs are able to apply the advantages of the central-
limit theorem in a variety of ways — they can produce regression adjusted
estimates (including incorporating doubly robust estimators), they can be
applied to maximum likelihood estimators, and even causal models.

Some experimentation solution providers (google_bayesian;
stucchio2015bayesian) have taken a different strategy, and have relied on
Bayesian Hypothesis tests (Chapter 11 of (wasserman2004all)). In addition, the
Bayes Factor has also been proposed as a way of enabling continuous monitoring
in A/B tests (deng2016continuous). Naturally, these represent radically
different philosophical approaches to analyzing A/B tests, and they provide a
different set of statistical guarantees. For the moment, we note that these
methods also rely on the specification of parametric models, and their
performance is generally sensitive to choices of priors (that technology
industry practitioners are generally not well equipped to determine).

###  2.1. Baseline Approaches & Evaluation

The hypothesis test of interest is the “difference of means” between the two
treatment arms (in later sections we extend to greater than two treatments).
The hypothesis being tested for an A/B test is (conducted at the level
$\alpha$):

(1) $H_{0}:\mu_{0}=\mu_{1}\quad\text{vs.}\quad H_{1}:\mu_{0}\neq\mu_{1}.$

The first method we consider is a classical z-test (in large sample limit,
this is equivalent to a t-test). Let $\hat{\theta}_{n}$ be the sample
difference in means, and $\hat{\sigma}_{n}$ be the sample (pooled) standard
deviation after $n$ samples have been recorded. For this test, the (point-
wise) confidence interval is given by
$CI_{1-\alpha}\coloneqq\hat{\theta}_{n}\pm
z_{1-\alpha/2}\frac{\hat{\sigma}_{n}}{\sqrt{n}}$. Note that the statement that
the null hypothesis is rejected (i.e., an effect size is non-zero) at
significance level $\alpha$ is equivalent to a statement that the $(1-\alpha)$
CI does not include the null. Formally, the type-I error guarantees of this
confidence interval are only satisfied if the result is looked at once, this
will be the first baseline approach and we will call it (1) FHT (fixed horizon
test without peeking).

In what follows, we will also explore an alternate stopping rule, the test
will be stopped as soon as the above $(1-\alpha)$ CI does not include the null
hypothesis. Of course, this is formally an incorrect way to run a test. But it
is emblematic of the peeking problem we are attempting to fix. This is our
second baseline, (2) FHT-peeking (fixed horizon test with peeking).

As our third baseline, we consider the group sequential method define by
(gordon1983discrete; demets1994interim), (3) LDM (Lan and De Mets). This
requires pre-specification of the total sample size, but allows flexible
peeking (with the requirement that the times of the peeks are pre-registered
and accounted for). We use the Pocock alpha-spending function
(pocock1977group), which spends $\alpha$ evenly over the length of the
experiment. The cumnulative $\alpha$ spent by the peek at time $t^{*}$ is
given by $\alpha\log(1+(e+1)t^{*})$, with the incremental $\alpha$ spent at
the $k^{\text{th}}$ peek given by $\alpha(t_{k}^{*})-\alpha(t_{k-1}^{*})$.

Our next baseline is the Mixture Sequential Probability Ratio Test ((4) mSPRT)
proposed in (johari2017peeking). A $(1-\alpha)$-CS that can be derived by
inverting the mSPRT is given by

$CS_{1-\alpha}^{\text{mSPRT}}\coloneqq\Biggl{\\{}\hat{\theta}_{n}\pm\sqrt{\frac{2(n\rho^{2}+\hat{\sigma}^{2}_{n})}{n^{2}\rho^{2}}\log\left(\sqrt{\frac{n\rho^{2}+\hat{\sigma}^{2}_{n}}{\alpha\hat{\sigma}_{n}}}\right)}\Biggr{\\}}.$

Here, $\rho^{2}$ is a free parameter that can be tuned. The test is stopped as
soon as the above $CS$ does not include the null hypothesis.

A decision theoretic way of running A/B tests is to instead define a loss
function $L(\theta,\theta^{*})$ for an unknown parameter $\theta$, which
quantifies the cost of deciding that the parameter has a different value
$\theta^{*}$ instead of its “true” value, $\theta$. While decision theory
spans both frequentist and Bayesian approaches, here we will focus on the
latter. Bayesian methods will generally give a posterior distribution
$p(\theta|X)$ for the unknown parameter $\theta$ after having observed data
$X$. Decision theory then has us compute the posterior expected loss
$E\left[L(\theta^{*}|X)\right]=\int L(\theta,\theta^{*})p(\theta|X)d\theta$,
which is a function of $\theta^{*}$. The optimal decision is then the choice
of $\theta^{*}$ that minimizes this posterior expected loss. For A/B testing,
there are of course many possible loss functions one can consider, e.g., the
company VWO, in a white paper on Bayesian methods (stucchio2015bayesian),
suggested a loss function that penalizes wrong choices linearly. For a binary
outcome, we implement (stucchio2015bayesian), with a uniform Beta prior for
propensity of response, so that the posterior for each (independent) arm is
also a Beta distribution. We evaluate this posterior expected loss at every
peek, and stop as soon as it is below a “threshold of caring,”
$\varepsilon$.This is our next baseline (5) BHT-uninformed.

Another popular paradigm for Bayesian “testing” is the use of Bayes Factors
(deng2016continuous). In this framework, we assume a prior probability for
$H_{i}$ to be true ($p(H_{i})$). After collecting data, $X$, for the
experiment, we compute the posterior odds of BF: $p(H_{1}|X)/p(H_{0}|X)$. The
BF is the ratio of the likelihood of observing data under $H_{1}$ relative to
$H_{0}$. If we assumed both the hypotheses to be equally plausible to start
with, the Prior Odds is $1$, and the Posterior Odds are equal to the Bayes
Factor. We then run the test, and stop it when the Posterior Odds are above
(or below) some threshold, e.g., above 10, (or below 0.1) representing 10:1
odds that the alternative (or null) is consistent with the data. In our
implementation, for binary outcomes, we assume that $H_{0}$ is a single Beta
distribution while for $H_{1}$, each arm is independently described by a Beta
distribution, all with uniform priors 111see details in supplementary material
online. We stop when the Posterior Odds are above $1/\alpha$ representing 20:1
odds that the alternative (or null) is consistent with the data given
$\alpha$=0.05. This is our final baseline, (6) BF-uninformed.

##  3\. Anytime Valid Confidence Sequences

###  3.1. Average Treatment Effect

The first consideration for an online A/B test is to define the appropriate
test statistic and the associated confidence bound for it. For this task, we
choose to consider the difference of treatment means (average treatment effect
or ATE) as the primary test statistic of interest to evaluate the hypothesis
test in Equation (1). That is, we are interested in the test statistic
($\hat{\theta}=\hat{\mu}_{1}-\hat{\mu}_{0}$).

The primary AsympCS defined in (waudby2021doubly) for the sample mean of a
random variable, requires independent and identically distributed observations
from a distribution with a finite $2+\delta$ moment ($\delta>0$).[We have
AsympCSs for data that are neither independent nor identically distributed.
–Ian] So, while we cannot directly apply their Theorem 2.2, we can use the
results derived for computing the CS for an ATE from observational causal
inference.[I’m confused by this — observational causal inference does not
imply non-iid. Indeed, the classical observational causal setup assumes iid
(even though we don’t). Also, why are we emphasizing _observational_ here? A/B
tests are randomized experiments, and thm 3.3 is exactly for those. –Ian]
Relying on Theorem 3.3 of this work, we model the influence functions of the
ATE, which relies of the inverse propensity weighted estimator (IPTW),
allowing us to write the ATE as the mean of functionals of the data. This
derivation leads to the following form of our confidence sequence with anytime
guarantees:

$\displaystyle
CS_{n}^{\text{Asymp}}\coloneqq\Biggl{\\{}(\hat{\mu}_{1;n}-\hat{\mu}_{0;n})\pm$
$\displaystyle\sqrt{\frac{n}{n-1}\left[\frac{n}{n_{0}}(\hat{\sigma}_{0;n}^{2}+\hat{\mu}_{0;n}^{2})+\frac{n}{n_{1}}(\hat{\sigma}_{1;n}^{2}+\hat{\mu}_{1;n}^{2})-(\hat{\mu}_{1;n}-\hat{\mu}_{0;n})^{2}\right]}\times$
(2)
$\displaystyle\sqrt{\frac{2(n\rho^{2}+1)}{n^{2}\rho^{2}}\log\left(\frac{\sqrt{n\rho^{2}+1}}{\alpha}\right)}\Biggr{\\}}.$

[I think it’s important for all of these estimators to have an $n$ in them to
emphasize that they change at each step. For example, $\widehat{\mu}_{n}^{1}$
instead of $\widehat{\mu}_{1}$ and so on. –Ian] Where, $n$ is the sample size
observed at the time of evaluation of the CS, of which $n_{0}$ and $n_{1}$ are
the samples on treatment $0$ and $1$, respectively. The running sample means
and standard deviations are given by $(\hat{\mu}_{0;n},\hat{\mu}_{1;n})$ and
$(\hat{\sigma}_{0;n},\hat{\sigma}_{1;n})$, respectively. The desired level of
the test is $\alpha$, and $\rho$ is a tuning parameter [A tuning parameter
that does what? Maybe we can just point the reader to Section C.3. –Ian]. The
above CS has the desirable property that $\mathbb{P}(\forall
n\in\mathbb{N}^{+},(\mu_{1}-\mu_{0})\in
CS_{1-\alpha}^{\text{Asymp}})>(1-\alpha)$ [not exactly true… –Aaditya]. [I
agree with Aadi, this is not true. This is an AsympCS, and hence it only
satisfies this guarantee in an approximate sense given by waudby2021doubly.
Also i would recommend removing $(1-\alpha)$ from the def of
$CS_{1-\alpha}^{\mathrm{Asymp}}$ and instead have a subscript $n$. Usually
$\alpha$ is chosen beforehand, so it’s ok to omit, but $n$ indexes the
different intervals that form the entire CS. –Ian]

###  3.2. Correction for Multiple Treatment Arms

Often, the experimenter is interested in comparing more than two treatment
arms (say $k$ arms). In such situations, we need to consider the multiple
hypothesis testing challenges that arise from having multiple treatments. We
avoid the possibility of testing all possible pairs from the $k$ arms, and
instead require the experimenter to define a control (or baseline) arm (say
$0$). Then for each of the $k-1$ other treatments we consider the hypothesis
$H_{0i}:\mu_{0}=\mu_{i}\text{ vs. }H_{1i}:\mu_{0}\neq\mu_{i}$, where
$i\in\\{1,2,\cdots,k-1\\}$. Applying multiple testing correction strategies
like Bonferroni (bonferroni1936teoria) or BH Procedure
(benjamini1995controlling) require the definition of $p$-values. While AsympCS
does not provide us with a $p$-value [We actually do have asymptotic
$p$-processes (asymptotic analogues of anytime $p$-values) defined in Def C.3.
Proposition C.3 gives an explicit one. This uses the duality between anytime
$p$-values and CSs. –Ian], we do have a definition of a confidence sequence.
Using the duality between a confidence bound and the $p$-value, we define the
following as our $p$-value.

Rewrite the confidence sequence in Equation (2) as
$CS_{1-\alpha}^{\text{Asymp}}\coloneqq\\{(\hat{\mu}_{1}-\hat{\mu}_{0})\pm\Gamma_{n}\\}$.
Then we arrive at the following

$\tilde{z}=\frac{\hat{\mu}_{1}-\hat{\mu}_{0}}{\Gamma_{n}/z_{1-\alpha/2}}\quad\text{and}\quad
p\text{-value}=2(1-\Phi(|\tilde{z}|)).$

[Where did this definition of a $p$-value come from? What does it satisfy,
etc.? I don’t think it’s making use of the duality. –Ian] Once we have the
$p$-values for the $k$ tests, we can directly apply multiple correction
procedures on them. For example, to apply the BH procedure, we do the
following. We rank order the $p$-values in increasing order, say
$q_{(1)},q_{(2)},\cdots,q_{(k-1)}$, and find the largest $j$ such that
$q_{(j)}<\alpha\cdot j/k$, and reject hypotheses $1,2,\cdots,j$. [Should
perhaps mention what guarantee this procedure satisfies. –Ian]

###  3.3. Test Statistic for Lift

Sometimes, an important marketing question is how much “lift” a treatment
generates over the control. Formally, this is defined as
$(\mu_{1}/\mu_{0}-1)$, and may also be represented as a percentage. Since this
is a commonly used test statistic in an A/B test, we would like to provide an
AsympCS for it. Unfortunately, since the natural estimate of this test
statistic $(\hat{\mu}_{1}/\hat{\mu}_{0}-1)$ cannot be written in the form of a
sum of random variables, we cannot directly apply the results of
(waudby2021doubly) to generate an AsympCS for this.

To generate an AsympCS for this test statistic, we do the following: (1) We
start the logarithm of the ratio of the treatment means. (2) We use the fact
that for a is $X$ is a continuous random variable, such that:
$\Pr(X<\ell)=\alpha/2$, then, for a strictly monotonic increasing function
$g:\mathbb{R}\rightarrow\mathbb{R}$, $\Pr(g(X)<g(\ell))=\alpha/2$. (3) Using
the one-sided AsympCS from (waudby2021doubly) and the monotonic transformation
$\log(\cdot)$, we generate a one-sided AsympCS for $\log\mu_{i}$. (4) to
generate the AsympCS for the difference $(\log\mu_{1}-\log\mu_{0})$, we take
the worst case bounds from each of the terms to arrive at the AsympCS for the
difference. (5) Finally, we apply another monotonic transformation to arrive
at the AsympCS for the lift test statistic.

The AsympCS for the lift can be shown to be

$CS_{1-\alpha}^{\text{Asymp-
Lift}}\coloneqq\Biggl{[}\frac{\hat{\mu}_{1}-\hat{\sigma}_{1}\beta(n_{1},\alpha,\rho)}{\hat{\mu}_{0}+\hat{\sigma}_{0}\beta(n_{0},\alpha,\rho)}-1,\frac{\hat{\mu}_{1}+\hat{\sigma}_{1}\beta(n_{1},\alpha,\rho)}{\hat{\mu}_{0}-\hat{\sigma}_{0}\beta(n_{0},\alpha,\rho)}-1\Biggr{]},$

where
$\beta(n_{0},\alpha,\rho)=\sqrt{\frac{2(n\rho^{2}+1)}{n^{2}\rho^{2}}\log\left(\frac{\sqrt{(n\rho^{2}+1)}}{\alpha}\right)}$,
and $(\hat{\mu}_{0},\hat{\mu}_{1})$ and $(\hat{\sigma}_{0},\hat{\sigma}_{1})$
are the running means and standard deviations in the two treatment arms. [This
is all really neat. I haven’t checked the derivations in detail but it looks
reasonable. Perhaps worth discussing at some other point: I think there’s a
better way to estimate lift, since it can be posed as a Z-estimation problem,
which has an influence function, and hence waudby2021doubly will apply. –Ian]

###  3.4. Sample Size Calculator

[Write sample size details –David]

##  4\. Experiments

###  4.1. Type I Error

First, we consider the type-I error of the baselines in Section 2.1 with our
proposed CS in Section 3.1. The considered type-I error is one where we
continuously monitor the test using the stopping rules from each method, when
there is no difference between the two treatments (i.e., the null hypothesis),
this may also be considered an “A/A test”, which is a common practice in
online experimentation to test the validity of the testing strategy. While not
all the baselines are designed to control this type-I error, given the
standard way in which A/B tests are analyzed in this industry, we strongly
feel that controlling this type-I error to the desired level of the test (say,
5%) is a highly valuable property.

Figure 2 presents the results from this analysis. All simulations are
performed for binary random variables. We see XXXX. [Update once the plot is
updated. –Ritwik]

Figure 2. PLACEHOLDE IMAGE, TO BE UPDATED - Cumulative miscoverage rates of
asymptotic confidence sequences of both mSPRT and AsympCS are bound at $0.05$.
However the credible intervals from the Bayesian approaches are inflated.

###  4.2. Statistical Power and Stopping Time

For the methods that maintain the desired Type-I error (LDM, mSPRT and
AsympCS), we next consider the statistical power of these tests. In Figure (3)
we plot the statistical power of these three methods against the multiple of
the FHT sample size that provides 80% power. We see that LDM ([How many peeks?
6 peeks –Simon]) has higher statistical power with the same sample size
compared to the anytime valid methods. We also note that AsympCS consistently
has higher power than mSPRT. When compared to the FHT sample size, we see that
to achieve an 80% power with AsympCS, we will need about two to three times
the sample size of an FHT with 80% power. For LDM, we need to store the times
of the ‘peeks’ and are unable to arbitrarily continue the test, given these
limitations, we favor the anytime valid methods in our implementation.

[I would like to remove the next parargraph - i dont think it’s a robustly
true statement –Akash] We also evaluate the stopping times of the test for
mSPRT and AsympCS, see Figure (4). We see that the stopping times for AsympCS
are consistently smaller than for mSPRT. We see that for the range of the
treatment mean, AsympCS takes $1.5$ to $4.3$ times of the FHT sample
size.222It is important to note that even though $\rho^{2}$ is a tuning
parameter in our approach, the results are robust to its choice. We perform
additional analysis, see Appendix, and choose a value that provides good
properties.

Figure 3. Statistical power for a binary response, as a function of the sample
size (multiple of the sample size for an FHT with 80% power). Figure 4.
Stopping time of the test as a function of the FHT sample size.

###  4.3. Sample Size Calculator

An analogue of fixed-horizon power can be derived by considering the expected
stopping time of an experiment. Specifically, given

###  4.4. Power of AsympCS for Lift

For the AsympCS for Lift, we evaluate the type I error and power. We observe
that this strategy has the correct type I error (not shown). We also explore
the power of using the proposed AsympCS for Lift. We see in Figure 5 that this
has lower power than the AsympCS for the difference of means. In spite of
that, we see that the power soon rises to $1$ as the true lift increases.

Figure 5. Power of AsympCS for the lift test statistic.

###  4.5. Real World Experiments

We consider $2,089$ real A/B tests conducted on our experimentation service.
As seen in our simulations, if we have relied on FHT with “peeking”, we would
allow the possibility erroneous conclusions and early stopping of the tests.
We observe a similar result in our real experiments. Based on the full length
of the conducted experiment, Table 1 shows that $57\%$ of the times both
strategies indicated that the treatment arms are not statistically significant
and $28\%$ of the times both agreed in concluding that the test is
statistically significant. Hence, there is agreement on the conclusion for
$85\%$ of the experiments. The most critical advantage by the implementation
of AsympCS is that $15\%$ of the time, it determined that the test is not
significant, while FHT with continuous monitoring inferred statistical
significance.

Table 1. Comparison of AsympCS and FHT on real A/B Tests

| AsympCS Significant | AsympCS Not Significant | Total
---|---|---|---
FHT Significant | 28% (593) | 15% (308) | 901
FHT Not Significant | 0.1% (3) | 57% (1185) | 1188
Total | 596 | 1493 | 2089

Figure 6. Sample size ratio.

A common criticism of anytime valid methods is that they take longer to reach
a conclusion compared to a correctly implemented FHT, i.e., single peek. But,
it is important to note that an FHT requires the specification of a “Minimum
Detectable Effect” (MDE). FHT with no peeks is efficient only if this MDE is
guessed correctly at the beginning of the test. If this MDE is underestimated,
the FHT is under-powered; if it is overestimated, then the test takes longer
to stop, and an anytime method could have reached the conclusion sooner. This
fact is displayed based on the real tests in Figure 6. We see that if the MDE
was correctly specified, we would in fact take longer to stop. But if the MDE
was 50% lower than specified, AsympCS would in fact stop with about about half
the samples for the FHT.333Additional details can be found in the Appendix.

###  4.6. Implementation in Adobe Experience Platform

Our solution is implemented in the Experimentation Service that is available
on Adobe’s Experience Platform. The AsympCS based analysis of online A/B tests
are available in the product Adobe Journey Optimizer (Getstart17:online;
Statisti60:online). In addition to that, the same technology is also available
in Customer Journey Analytics (cja), which is a flexible data analysis
ecosystem which customers can connect to their own experimental data and
analyse with our proposed approach. [To update –Akash]

##  5\. Conclusion

We have proposed and implemented an anytime-valid A/B testing strategy for
online A/B tests based on asymptotic confidence sequences, permitting
continuous monitoring of experimental results. Our approach has anytime
guarantees, and performs better than other anytime-valid methods, in terms of
statistical power. Our implementation, in Adobe’s Experience Platform,
additionally allows handling multiple treatments, alternate test statistics
like lift, and provides a handy sample size calculator to design A/B tests. We
show the advantages of our approach using simulations, and on several real A/B
tests conducted on our platform. In the future, it will be interesting to
explore anytime valid methods for multi-arm bandits.


# Minor Issues Escalated to Critical Levels in Large Samples: A Permutation-
Based Fix

Xuekui Zhang
Department of Mathematics and Statistics, University of Victoria, Canada
and
Li Xing
Department of Mathematics and Statistics, University of Saskatchewan, Canada
and
Jing Zhang
Department of Mathematics and Statistics, University of Victoria, Canada
and
Soojeong Kim
Department of Mathematics and Statistics, University of Victoria, Canada This
research was funded by NSERC DG: # RGPIN-2022-03050 (XZ), # RGPIN-2021-03530
(LX) and the Canada Research Chair #CRC-2021-00232 (XZ), Michael Smith Health
Research BC Scholar: # SCH-2022-2553(XZ). This research was enabled in part by
computational resource support provided by Westgrid (https://www.westgrid.ca)
and the Digital Research Alliance of Canada (https://alliancecan.ca).

###### Abstract

In the big data era, the need to reevaluate traditional statistical methods is
paramount due to the challenges posed by vast datasets. While larger samples
theoretically enhance accuracy and hypothesis testing power without increasing
false positives, practical concerns about inflated Type-I errors persist. The
prevalent belief is that larger samples can uncover subtle effects,
necessitating dual consideration of p-value and effect size. Yet, the
reliability of p-values from large samples remains debated.
This paper warns that larger samples can exacerbate minor issues into
significant errors, leading to false conclusions. Through our simulation
study, we demonstrate how growing sample sizes amplify issues arising from two
commonly encountered violations of model assumptions in real-world data and
lead to incorrect decisions. This underscores the need for vigilant analytical
approaches in the era of big data. In response, we introduce a permutation-
based test to counterbalance the effects of sample size and assumption
discrepancies by neutralizing them between actual and permuted data. We
demonstrate that this approach effectively stabilizes nominal Type I error
rates across various sample sizes, thereby ensuring robust statistical
inferences even amidst breached conventional assumptions in big data.
For reproducibility, our R code will be publicly available at:
https://github.com/ubcxzhang/bigDataIssue before April 2024.

Keywords: big data, hypothesis testing, inflated Type I error, violated model
assumptions

##  1 Background

In the era of big data, analyzing datasets with vast volumes of entries
prompts a reevaluation of the applicability of traditional statistical methods
to such extensive data. As a basic fact stated in statistics textbooks, larger
sample sizes are theoretically favoured for increased accuracy, reduced
variability, and enhanced hypothesis testing power without intensifying Type I
error rates. However, in practical scenarios, especially in medical studies,
there’s a concern that large sample sizes could inflate Type I errors, a point
not fully explored regarding underlying causes (Armstrong, 2019; Faber and
Fonseca, 2014). A widely accepted view is that larger samples can reveal even
minor effects, making it crucial to consider both p-value and effect size in
large datasets (Dunkler et al., 2020; Sullivan and Feinn, 2012). However, the
debate (on forums like _StackExchange_) continues over whether huge samples in
practice inflate traditional methods’ p-values, which is the focus of this
work.

Our experiences in differential expression (DE) analysis of single-cell RNA-
sequencing (scRNA-seq) data underscore the prevalent concerns regarding
inflated p-values. DE analysis, a fundamental aspect of genomic data analysis,
aims to identify genes with significantly different expression levels under
various conditions (McDermaid et al., 2018). The scRNA-seq is a high-
throughput technique capturing gene expression patterns at the granularity of
individual cells. Over the past eight years, scRNA-seq sample sizes have
escalated dramatically, growing from a single cell to hundreds of thousands
(Svensson et al., 2018). Most DE methods developed for scRNA-seq data emerged
at the inception of this new technology when the number of cells was either
small or moderate. Yet, contemporary scRNA-seq data analysis with these DE
methods identifies thousands of significant genes (even after multiple testing
adjustments), resulting from a large sample size. This has led researchers to
rely on subjective log-fold change thresholds (such as 0.5, 1, or 2) to
discern differentially expressed genes, presenting two critical problems.
First, given the relationship between effect sizes and p-values, the p-values
are nearly obsolete for identifying DE genes. Second, it raises the question
of whether these numerous small p-values indicate minor effects or if some are
misleadingly insignificant due to the large sample sizes involved.

We argue that large samples don’t inherently cause issues, but can amplify
pre-existing minor problems to critical levels, demanding careful handling. In
practice, the ideal conditions posited by statistical theory are rarely
encountered, leading to complex issues. Such problems, stemming from unmet
model assumptions or data collection flaws, may become more pronounced with
larger samples, indicating the necessity for novel strategies to tackle these
challenges. To address these challenges, we propose a permutation-based method
for deriving p-values unaffected by sample size, even in non-ideal conditions.
We utilize a permutation approach since it is flexible to be applied to
samples from any distributions (Pitman, 1937, 1938).

Our simulation study illustrates how sample size affects Type I errors in
conventional methods and demonstrates our approach’s corrective capabilities.
In this study, we deliberately violate model assumptions in two common real-
world scenarios. Firstly, we create data with the outcome variable following a
distribution different from the model assumption, inspired by the fact that
all distribution assumptions in DE analysis over-simplify the intricate nature
of real-world data. Secondly, we simulate outcome variable values using two
predictors but only include one in the analysis, reflecting situations where
not all influencing factors are measurable or considered. This approach
mirrors the complexities often encountered in real-world scenarios where
comprehensive data capture is challenging.

##  2 Notations and Problem Setup

To investigate how large sample size affects Type I error of statistical
analysis with a simulation study, we consider the problem of fitting a Poisson
regression to analyze a dataset comprising one outcome variable and one
predictor. The objective is to determine the statistical significance of the
predictor’s association with the outcome variable, primarily through the
p-value of the regression coefficient for the predictor. The model is
structured as follows:

$\displaystyle y_{i}$
$\displaystyle\sim\mathrm{Poisson}(\lambda_{i}=e^{\beta_{0}+\beta_{1}x_{1i}})$
(1)

Here, $i=1,\ldots,n$ denotes the sample index. The term $y_{i}$ represents the
outcome following a Poisson distribution with parameter $\lambda_{i}$, and
$x_{1i}$ is the predictor for the $i$-th sample.

In situations involving large sample sizes, the concern for statistical power
diminishes, shifting our attention primarily to the Type I error rate. To
investigate this, we set $\beta_{1}=0$ and focus on the proportion of false
positives detected in our analysis. Our simulation study considers two
specific scenarios, representing two common issues encountered in real-world
data analysis.

Firstly, we explore a scenario of distribution misspecification. Here, rather
than using a Poisson distribution, the outcomes $y_{i}$ are derived from a
discretized $F$ distribution with degrees of freedom $(8,8)$. These values are
subsequently rounded to the nearest integer:

$\displaystyle z_{i}$ $\displaystyle\sim F_{8,8}\quad\mathrm{and}\quad
y_{i}=\mathrm{round}(z_{i})$ (2)

Secondly, we consider a scenario involving unobserved predictors. In this
case, the outcome $y_{i}$ is influenced by an unobserved predictor $x_{2i}$,
modelled as follows:

$\displaystyle y_{i}$
$\displaystyle\sim\mathrm{Poisson}(\lambda_{i}=e^{\beta_{0}+\beta_{2}x_{2i}}).$
(3)

For each scenario, various settings with differing parameter values will be
examined. To estimate the Type I error for each setting, we simulate $K=1000$
datasets and fit the Poisson regression model as per (1). Denote the p-values
for $\beta_{1}$ as $p_{j}$ for $j=1,\ldots,K$. The Type I error rate for a
given setting is then defined as:

$\displaystyle\mathrm{Type\ I\ error\ rate}$
$\displaystyle=\frac{\sum_{j=1}^{K}I(p_{j}<0.05)}{K}$ (4)

When a method accurately maintains the Type I error rate at the nominal level
of 5%, the decision to reject the null hypothesis in any given simulated
dataset follows a Bernoulli distribution with the parameter $p=0.05$.
Consequently, the total count of rejections across $K=1000$ repeated
experiments adheres to a binomial distribution, specifically
Binom$(1000,p=0.05)$. This binomial distribution yields a 95% confidence
interval for the estimated Type I error, ranging between (0.0362, 0.0638),
which provides a statistical measure of the variability expected around the
nominal 5% Type I error rate. Our simulation study aims to assess the impact
of increasing sample sizes on the defined Type I error.

##  3 Simulation Study

###  3.1 Impact of Large Sample Size on the Estimated Parameters of a
Distribution

Before tackling the Poisson regression challenge, we first investigate a
simpler problem to demonstrate how increased sample sizes can influence the
parameter estimation of a Poisson distribution when data collection is flawed.
We consider a Poisson distribution where each event count $y_{i}$ follows
$y_{i}\sim\mathrm{Poisson}(\lambda=5)$ for $i=1,\ldots,n$. However, due to the
limited sensitivity of the measuring device, any counts less than 2 are
erroneously recorded as 0. This situation sets the stage for understanding the
repercussions of discrepancies between the assumed and actual distributions as
the number of observations, $n$, increases.

We anticipate that the measuring device’s insensitivity will introduce a bias
in estimating $\lambda$, and we aim to analyze how this bias shifts with
enlarging sample sizes. To do this, we simulate datasets of 10 varying sample
sizes equally spaced in log-10 scale, i.e. $\mathrm{log}_{10}(n)$ takes 10
values equally spaced between 1 and 6. Each dataset is generated under the
flawed measurement condition as described above; from these, we estimate the
Poisson parameter $\hat{\lambda}=\sum_{i=1}^{n}y_{i}/n$. We focus on examining
the relationship between the sample size $n$ and the estimation bias
$\hat{\lambda}-\lambda$. To assess the uncertainty in our estimates, we
replicate the simulation 1000 times for each sample size, yielding 1000 bias
measurements per size.

Figure 1: Each boxplot displays the distribution of 1000 biases linked to the
estimated Poisson rate parameter $\lambda$ for a specific sample size,
indicated on the y-axis. The x-axis denotes the magnitude of these biases,
with a dashed red line at 0 indicating the point of unbiased estimation. The
boxplots reveal a decrease in estimation uncertainty as the sample size grows,
indicated by the shrinking interquartile ranges (IQRs) or the size of boxes.
Nevertheless, while precision improves, it also introduces a systematic
positive bias that consistently emerges with larger sample sizes.

Each boxplot in Figure 1 displays the distribution of 1000 biases associated
with the estimated Poisson rate parameter $\lambda$ for a specific sample
size. The x-axis denotes the magnitude of these biases, with a dashed red line
at 0 indicating the point of unbiased estimation. The boxplots reveal a
decrease in estimation uncertainty as the sample size grows, indicated by the
shrinking interquartile ranges (IQRs) or the size of boxes. However, this
increased precision comes with a systematic positive bias that persists across
larger sample sizes. In more complicated analyses, this bias in parameter
estimation could be consequential, affecting the accuracy of confidence
intervals and p-values. For smaller sample sizes, the wider confidence
intervals tend to mitigate the effects of bias, reducing its influence on the
analysis. However, as sample sizes increase and confidence intervals
consequently narrow, even small systematic biases can become critical, leading
to substantial errors in inferential decision-making. This emphasizes the
importance of accounting for and addressing bias, especially in large-scale
data analyses. Next, we will show how growing the sample size affects Poisson
regression.

###  3.2 Impact of Large Sample Size on Type I Error in Poisson Regression
with Assumption Violations

We now return to the primary focus of this simulation study: Poisson
regression. Hypothesis testing of Poisson regression’s coefficient is more
complex than the basic estimation of a Poisson distribution’s parameter, as it
involves multiple variables (especially with more than one predictor) and
necessitates the determination of a null distribution for the test statistics
of the regression coefficients.

Using simulation studies, we explore two distinct scenarios of assumption
violations previously outlined in the model framework section. Here, we fit a
Poisson regression model, as defined in Equation (1), but model assumptions
are violated in two different ways. For both scenarios, $y$ and $x_{1}$ are
generated as unrelated by setting $\beta_{1}=0$, making any significant
p-values indicative of Type I errors instead of detecting small effects. We
aim to examine the influence of increasing sample size on Type I error by
analyzing various settings with 60 different sample sizes. Specifically, we
let $\mathrm{log}_{10}(n)$ takes 30 values equally spaced between 1 and 2, and
another 30 equally spaced between 2 and 5. For each sample size setting, we
simulate $K=1000$ datasets to estimate the corresponding Type I error rate
using the formula provided in Equation (4). More densed sample sizes are
investigated in the range $(10^{1},10^{2})$ than the range $(10^{2},10^{5})$
for two reasons: (1) smaller sample sizes take less computing resources and
(2) type I error rate charge faster in the first region in our simulation
study.

Figure 2: The relationship between the sample sizes (x-axis) and the Type I
error rates (y-axis), demonstrating how the issue of misspecified outcome
distribution is amplified by growing sample sizes. Each circle represents a
Type I error rate of Poison regression estimated from $K=1000$ simulated
datasets, and the solid curve represents smoothed pattern of these circles.
Each cross represents a Type I error rate of permutation-corrected Poison
regression estimated from the same $K=1000$ simulated datasets, and the dashed
curve represents smoothed pattern of these crosses. Shaded area represents the
95% confidence interval of observed Type I error rate if the method can
control it at the nominal level.

In the first scenario, the actual distribution of the outcomes diverges from
the Poisson distribution that the model presumes. Specifically, the outcomes
$y_{i}$ are generated from a discretized $F$ distribution with degrees of
freedom $(8,8)$, as detailed in Equation (2). Predictor data $x_{1i}$ are
drawn from a standard normal distribution $N(0,1)$. Figure 2 illustrates the
relationship between the sample sizes (x-axis) and the Type I error rates
(y-axis), demonstrating how the issue of misspecified outcome distribution is
amplified by growing sample sizes. Each circle represents a Type I error rate
of Poison regression estimated from $K=1000$ simulated datasets, and the solid
curve represents smoothed pattern of these circles. Shaded area represents the
95% confidence interval of observed Type I error rate if the method can
control it at the nominal level. The estimated Type I error rates of Poison
regression remain within the 95% confidence bounds when the sample size is
small. As sample size increases Type I error rates start to fluctuate in and
out of the confidence bounds, reflecting inconsistent control over false
positives. When the sample size keep growing, the estimated Type I error
predominantly lies outside these bounds, indicating that larger samples
exacerbate Type I error inflation when models use wrong distribution
assumptions. The Type I error rates become stable at certain values when
sample size is large enough, which converges to the theoratical value depends
on how the model assumption is voilated.

Both the discretized $F$ distribution and the Poisson distribution are
integer-valued and skew towards the right, possessing heavier tails on that
side. Our simulation study indicates that when the sample size is small, a
misspecified distribution poses a minor problem, and approximating a Poisson
distribution with a discretized F distribution is generally acceptable.
However, as the sample size increases, the ability to differentiate between
the two distributions enhances, thereby aggravating the issue of distribution
misspecification. This escalation becomes critical, leading to incorrect
decisions, notably inflated Type I errors.

In the second scenario, outcomes influenced by an unobserved predictor.
Specifically, $y_{i}$ are generated using $x_{2}$ according to the true model
specified in Equation (3). Consequently, $y$ and $x_{1}$ share no
relationship, and the actual predictor $x_{2}$ is not included in the data
analysis model described in Equation (1). In this setup, we simulate both the
noise predictor $x_{1i}$ and the genuine predictor $x_{2i}$ from a standard
normal distribution $N(0,1)$. We explore four different settings of true
regression coefficient values in model (3), particularly selecting
combinations of $\beta_{0}=0.3,0.5$ controlling the magnitude of baseline
Poison rates and $\beta_{2}=0.7,0.8$ controlling the strength of unobserved
signal. Figure 3 illustrates the relationship between the sample sizes
(x-axis) and the Type I error rates (y-axis). Same as in the first scenario,
each circle represents a Type I error rate of Poison regression estimated from
$K=1000$ simulated datasets, and the solid curve represents smoothed pattern
of these circles. Shaded area represents the 95% confidence interval of
observed Type I error rate if the method can control it at the nominal level.
Four coefficient settings are represented by different colors, as annotated in
the figure legend. Notably, we observed that the Type I errors in all settings
exceeded the 95% confidence interval’s upper bound 0.0638, indicating a severe
inflation. In this scenario, the distribution of outcomes $y_{i}$ is
influenced by two factors: the actual counts that follow a Poisson
distribution and Gaussian noise attributable to the unobserved $x_{2}$, which
is not accounted for in the model and consequently contributes to the residual
variation. This omission of a variance source leads to overconfident
inferences, which explains the exceedingly high Type I errors, even with as
few as 10 samples. The simulation study reveals that the impact of this
omission intensifies as the sample size increases up to $n=100$. Beyond this
point, the impact begins to plateau, ultimately aligning with the theoretical
implications of disregarding a normal noise component in the analysis. Upon
comparing results across four parameter settings, we observed that higher
values of $\beta_{2}$—indicative of stronger missing signals—correspond to
more severe inflation of Type I error rates. This is consistent with the
notion that the greater the unobserved information, the more compromised the
analysis results are. The variance in the data originates from two components:
the inherent randomness in the Poisson distribution, governed by $\beta_{0}$,
and the variation in the Normal distribution of unobserved predictor, dictated
by $\beta_{2}$. Poisson regression uses Poison distribution to encapsulate
uncertainty and fails to distinguish between these two noise sources.
Consequently, we also observed that increased $\beta_{0}$ values, which
amplify the variance attributed to the Poisson component, further exacerbate
the inflation of Type I error rates when there is an unobserved predictor.

Comparing results among four parameter settings, we found larger $\beta_{2}$
values (i.e. missing stronger signals) lead to more severe inflation of Type I
error rates, which agrees with our intuition that the more unobserved
information, the worse analysis results. The variance in data comes from two
sources: the randomness in Poison distribution controlled by parameter
$\beta_{0}$ and in Normal distribution controlled by $\beta_{2}$. When the
model only uses Poisson distribution to model the uncertainty, noises from two
sources cannot be isolated; hence, we also found that larger $\beta_{0}$
values (increasing Poison variance) lead to more severe inflation of Type I
error rates.

Figure 3: The relationship between the sample sizes (x-axis) and the Type I
error rates (y-axis), demostrating how the impact of missing a predictor can
be amplified by growing sample sizes. Each circle represents a Type I error
rate of Poison regression estimated from $K=1000$ simulated datasets, and the
solid curve represents smoothed pattern of these circles. Each cross
represents a Type I error rate of permutation-corrected Poison regression
estimated from the same $K=1000$ simulated datasets, and the dashed curve
represents smoothed pattern of these crosses. Shaded area represents the 95%
confidence interval of observed Type I error rate if the method can control it
at the nominal level. Four coefficient settings are represented by different
colors, as annotated in the figure legend.

##  4 Using Permutation-Based Methods to Address Issues Exacerbated by Large
Sample Sizes

Our simulation studies demonstrate that large sample sizes can inflate the
Type I error rate when the model assumptions of Poisson regression are
violated, either due to outcome distribution misspecification or the presence
of unobserved predictors. When the sample size is large enough, the damage of
assumption violations converges to its theoretical level. These issues are
rooted in an incorrect theoretical null distribution. To counteract these
problems, we propose a permutation-based approach to construct a robust null
distribution for the test statistic, specifically $\hat{\beta}_{1}$, the
estimate from Poisson regression. This method aims to ascertain the
distribution of $\hat{\beta}_{1}$ under the null hypothesis of no relationship
between $x_{1}$ and $y$.

Initially, we fit a Poisson regression model to the original data to obtain
the observed coefficient, $\hat{\beta}_{1}^{orig}$. We then create numerous
permuted versions of the original dataset (e.g., $N=1000$) by randomly
shuffling the values of $\mathbf{x}_{1}$ to break any actual connection
between $x_{1}$ and $y$. For each permutation, we re-estimate
$\hat{\beta}_{1}$ after fitting the Poisson regression model to the permuted
data, generating a new coefficient, $\hat{\beta}_{1,perm}$. These coefficients
from all permutations form an empirical null distribution for
$\hat{\beta}_{1}$ under the hypothesis of no association.

Upon completing all permutations, the p-value is computed as the proportion of
permuted datasets where the absolute value of the permuted coefficient equals
or exceeds that of the original coefficient. This p-value indicates the
probability of observing the given association between $x_{1}$ and $y$ if, in
fact, no such association exists.

This permutation-based method offers robustness against model misspecification
and unobserved predictors. Issues affecting the original dataset equally
impact the permuted datasets used to construct the null distribution.
Moreover, as the sample size increases, it does so for both the original and
permuted data, thereby neutralizing the effects of assumption violations and
sample size influences. The entire process is succinctly encapsulated in
Algorithm 1, offering a systematic approach to addressing the challenges of
large sample sizes in statistical analysis.

Algorithm 1 Permutation-based p-value for Poisson Regression

1:$\mathbf{y}=(y_{1},\ldots,y_{n})$, $\mathbf{x}_{1}=(x_{1i},\ldots,x_{1n})$,
$N$ (default $N=1000$)

2:Fit the Poisson regression model
$y\sim\mathrm{Poisson}(\lambda=e^{\beta_{0}+\beta_{1}x_{1}})$ to the data
$\mathbf{y}$ and $\mathbf{x}_{1}$.

3:Record the estimate $\hat{\beta}_{1}^{orig}$ from the original data.

4:Initialize a counter $count=0$ for the number of times permuted $\beta_{1}$
exceeds the original.

5:for $j=1$ to $N$ do

6: Permute $x_{1i}$ to get $x_{1,perm}^{(j)}$.

7: Fit the model to the permuted data
$y\sim\mathrm{Poisson}(\lambda=e^{\beta_{0}+\beta_{1,perm}x_{1,perm}^{(j)}})$.

8: Compute the estimate $\hat{\beta}_{1,perm}^{(j)}$ from the permuted data.

9: if $|\hat{\beta}_{1,perm}^{(j)}|\geq|\hat{\beta}_{1}^{orig}|$ then
$count\leftarrow count+1$.

10:Compute the p-value as $p=count/N$.

11:return p-value $p$.

We implemented the permutation-based algorithm on data generated in our
simulation study. The crosses and dashed curves in Figures 2 and 3 show
estimated Type I errors when using our permutation-based approach to analyze
the same datasets generated in the Simulation Section. We discovered that our
method effectively corrects the inflated Type I errors of the original Poisson
regression to mostly within the shaded area representing 95% confidence
interval the nominal level of $0.05$. Furthermore, the dashed trend curves of
all four settings fluctuate around 0.05 with no discernible trend between the
sample sizes and the Type I errors rates obtained through permutation
correction. This indicates the robustness of this method in maintaining error
rates in big data.

##  5 Conclusion and Discussion

In our simulation studies, we opted to permute $x_{1}$ rather than $y$, a
decision that is not critical when the analysis involves only two variables.
However, when multiple predictors are incorporated into the regression
analysis, permuting $x_{1}$ becomes essential. This is because our primary
objective is to scrutinize the relationship between $x_{1}$ and $y$, while
accounting for the effects of other predictors. Permuting $y$ in such contexts
would disrupt its relationship not only with $x_{1}$ but also with the other
covariates. Consequently, the null distribution derived from permuting $y$
would no longer accurately reflect the true relationships between $y$ and the
other covariates when it’s assumed to be unrelated to $x_{1}$. Therefore, the
specific choice of permutation is crucial in maintaining the integrity and
interpretability of the analysis, especially in more complex regression models
with multiple predictors.

In our simulation studies, we observed that the detrimental effects of model
assumption violations become more pronounced with increasing sample sizes, yet
the extent of these damages diminishes as the sample size continues to grow.
This observation parallels the well-known phenomenon where the benefits of
enhanced statistical power from increasing sample size diminish at larger
scales. Both phenomena stem from the same underlying principle: as the sample
size approaches infinity, the information extracted from the data
asymptotically approaches the theoretical truth of the data’s characteristics.
Consequently, both the benefits and the detriments of analysis are bound to
converge towards their theoretical limits. Therefore, when the current value
is near this limit, further enlarging the sample size yields minimal changes.
This principle underscores the diminishing returns of increasing sample size,
both in terms of benefits and in terms of mitigating the effects of assumption
violations.

The Poisson regression model used in our simulation study exemplifies how
certain statistical analyses necessitate specialized methods to uphold nominal
Type I error rates when dealing with large-sample data. Nevertheless, it’s
important to recognize that not all analytical methods require such
permutation-based corrections. For instance, the Central Limit Theorem (CLT)
posits that the sample mean converges to a normal distribution as the sample
size increases, assuming certain conditions are met (Cam, 1986). Consequently,
if the null distribution of a test is derived based on the assumption that the
sample mean follows a normal distribution, such methods do not require
correction in large-sample data analysis. This highlights the nuanced nature
of statistical methodologies, where the need for adjustments like permutation-
based remedies is contingent upon the specific assumptions and characteristics
of the test in question.

In conclusion, the correct p-value from some statistical methods relies on
correctly specifying the model and including all relevant variables. Sample
size increase may exacerbate the problem caused by incorrect model assumptions
or flaws in data collection. So, when analyzing large sample data, special
treatment may be necessary to obtain the correct p-values. We proposed a
permutation-based p-value that is robust to many types of model
misspecification and unobserved confounders, and its p-values are unaffected
by sample sizes. This makes it a powerful tool, especially in complex real-
world scenarios where assumptions about the data or model may not fully hold
and the sample size is large.

## References

  * Armstrong (2019) Armstrong, R. A. (2019).  Is there a large sample size problem?  Ophthalmic and Physiological Optics 39(3), 129–130.
  * Cam (1986) Cam, L. L. (1986).  The Central Limit Theorem Around 1935.  Statistical Science 1(1), 78–91.
  * Dunkler et al. (2020) Dunkler, D., M. Haller, R. Oberbauer, and G. Heinze (2020).  To test or to estimate? P‐values versus effect sizes.  Transplant International 33(1), 50–55.
  * Faber and Fonseca (2014) Faber, J. and L. M. Fonseca (2014).  How sample size influences research outcomes.  Dental Press Journal of Orthodontics 19(4), 27–29.
  * McDermaid et al. (2018) McDermaid, A., B. Monier, J. Zhao, B. Liu, and Q. Ma (2018).  Interpretation of differential gene expression results of RNA-seq data: review and integration.  Briefings in Bioinformatics 20(6), 2044–2054.
  * Pitman (1937) Pitman, E. J. G. (1937).  Significance Tests Which May be Applied to Samples from Any Populations. II. The Correlation Coefficient Test.  Supplement to the Journal of the Royal Statistical Society 4(2), 225–232.
  * Pitman (1938) Pitman, E. J. G. (1938).  Significance Tests which May be Applied to Samples from any Populations: III. The Analysis of Variance Test.  Biometrika 29(3/4), 322.
  * Sullivan and Feinn (2012) Sullivan, G. M. and R. Feinn (2012).  Using Effect Size—or Why the P Value Is Not Enough.  Journal of Graduate Medical Education 4(3), 279–282.
  * Svensson et al. (2018) Svensson, V., R. Vento-Tormo, and S. A. Teichmann (2018).  Exponential scaling of single-cell RNA-seq in the past decade.  Nature Protocols 13(4), 599–604.