For each hypothesis testing setting where one would normally use a p-value, a savi test can be designed, with a number of advantages that are elaborately described and illustrated in this vignette. Currently, this package includes savi tests for the z-test, t-test, Fisher’s exact test, the chi-squared test (the savi test of two proportions), and the logrank test for survival data. In this vignette, we will illustrate the concepts of savi testing and e-values with the t-test as a running example. This savi test is designed to, on average, detect the effect as quickly as possible, if the effect actually exists.

A big advantage of e-values over their p-value equivalents is that savi tests conserve the type I error guarantee (false positive rate) regardless of the sample size. This implies that the evidence can be monitored as the observations come in, and the researcher is allowed to stop the experiment early (optional stopping) without over-inflating the chance of a false discovery. By stopping early fewer participants will be put at risk. In particular, those patients who are assigned to the control condition, when a treatment is effective. Savi tests also allow for optional continuation, that is the extension of an experiment regardless of the motivation. For instance, if more funds become available, or if the evidence looks promising and the funding agency, a reviewer, or an editor urges the experimenter to collect more data.

Importantly, for the savi tests presented here neither optional stopping nor continuation leads to the test exceeding the tolerable type I error $\alpha$. Savi tests allow for anytime valid inferences, because the results do not depend on the planned, current, or future sample sizes. We illustrate these properties below.

Firstly, we show how to design an experiment based on savi tests for testing means.

Secondly, simulations are run to show that savi tests indeed conserve the type I error guarantee under optional stopping for testing means. We also show that optional stopping causes the false null rejection rate of the classical p-value test to exceed the tolerable level $\alpha$ type I error guarantee. In other words, with classical tests one cannot adapt to the information acquired during the study without increasing the risk of making a false discovery.

Lastly, it is shown that optionally continuing non-significant experiments also causes the p-value tests to exceed the tolerable level $\alpha$ type I error guarantee, whereas this is not the case for savi tests.

This demonstration further emphasises the rigidity of experimental designs when inference is based on a classical test: the experiment cannot be stopped early, nor extended. Thus, the planned sample size has to be final. As such, a rigorous protocol needs to account for possible future sample sizes, which is practically impossible. Even if such a protocol can be made, there is no guarantee that the experiments go exactly according to plan, as things might go wrong during the study.

The ability to act on information that accumulates during the study –without sacrificing the correctness of the resulting inference– was the main motivation for the development of savi tests, as it provides experimenters with the much needed flexibility.

Installation

The stable version can be installed by entering in R:

install.packages("safestats")

The development version can be found on GitHub, which can be installed with the remotes package from CRAN by entering in R:

remotes::install_github("AlexanderLyNL/safestats", build_vignettes = TRUE)

The command

library(safestats)

loads the package. We use the following colours in our plots:

freqColours <- c("#E31A1CE6", "#FB9A9980")
eColours <- c("#1F78B4E6", "#A6CEE380")
eColoursAlt <- c("#B15928E6", "#FFFF9980")

1. Designing savi t-test experiments

Type I error and power

To avoid bringing an ineffective medicine to the market, experiments need to be conducted in which the null hypothesis of no effect is tested. Here we show how flexible experiments based on savi tests can be designed.

As the problem is statistical in nature, due to variability between patients, we cannot guarantee that all of the medicine that pass the test will indeed be effective. Instead, the target is to bound the type I error rate by a tolerable $\alpha$ level, typically, $\alpha = 0.05$. In other words, at most 5 out of the 100 ineffective drugs are allowed to pass the test.

At the same time, we would like to detect an effect with high chance, say, 80% power, if it is present.

Case (I): Designing experiments where the minimal clinically relevant effect size and the targeted power are known. Target: The sample size to plan for

Not all effect sizes are equally important, especially, when a minimal clinically relevant effect size can be formulated. For instance, suppose that a population of interest has a population average systolic blood pressure of $\mu = 120$ mmHg (milimetre of mercury) and that the population standard deviation is $\sigma = 15$. Suppose further that all approved blood pressure drugs change the blood pressure by at least 9 mmHg, then a minimal clinically relevant mean difference is $\psi_{\min}=\mu_{\text{pre}} - \mu_{\text{post}} = 12$, and the minimal clinically relevant standard effect size is $\delta_{\min} = (\mu_{\text{pre}} - \mu_{\text{post}}) / (\sqrt{2} \sigma) = 12 / (15 \sqrt{2} ) = 0.57$, where $\mu_{\text{pre}}$ represents the average blood pressure before treatment and $\mu_{\text{post}}$ the average blood pressure after treatment of the population of interest.

With an acceptable type I error rate set at $\alpha=0.05$, a desired power of 80%, and a minimal clinically relevant effect size of $\delta_{\min}=0.57$, our objective is to design an experiment with a planned sample size that is able to detect $\delta_{\min}=0.57$. As we will see below, this planned sample size (nPlan) is non-strict, which makes the experiment flexible. To obtain the planned sample size we run the following code:

alpha <- 0.05
power <- 0.8
deltaMin <- 12/(sqrt(2)*15)
sigma <- 15

designObj <- designSaviT(deltaMin=deltaMin, alpha=alpha,
                         power=power, sigma=sigma,
                         alternative="greater", 
                         testType="paired", seed=1, pb=FALSE)

designObj
#> 
#>  Savi Paired Sample T-Test Design
#> 
#>               n1Plan±2se, n2Plan±2se = 36±2.147258, 36±2.147258
#>               n1Mean±2se, n2Mean±2se = 22±0.6510288, 22±0.6510288
#> minimal standardised mean difference = 0.5656854
#>                          alternative = greater
#>                                power = 0.8
#>                      e-variable type = mom
#>                      parameter: gMom = 0.16
#>                                alpha = 0.05
#>       Reject H_0: e-value >= 1/alpha = 20
#> 
#> Timestamp: 2026-09-05 21:09:00 CEST
#> 
#> Note: If it is only possible to look at the data once, then n1Plan = 44 and n2Plan = 44.

The design object defines both the parameter gMom that will be used to compute the e-value, e.g. 0.16, and the planned sample size(s) under optional stopping, e.g. 36, 36. Hence, in this case we need the pre- and post-measurements of about 36 patients to detect a true standarised effect size of $\delta_{\min} = 0.57$. This nPlan of 36 is based on continuously monitoring the e-value and stopping the experiment as soon as it exceeds $1/\alpha = 20$. The first time/sample size that the E-variable hits the threshold $1/\alpha = 20$ is data dependent, thus, random. This randomness is expressed with nPlan being reported with two (bootstrap) standard error of the mean. When it is only possible to conduct the test once, when the data are treated as a single batch, then 44 patients (thus 8 more) are needed to detect $_{} = 0.57 $ with 80% chance.

Case (II): Minimal clinically relevant effect size and a budget constrained number of samples known. Target: Power

The following scenario involves a known minimal clinically relevant effect size, and a budget constraint such that we can at most invite, say, 30 participants for our study. The power of the test can then be explored with the following code:

designObj2 <- designSaviT(deltaMin=deltaMin, alpha=alpha,
                          nPlan=c(30, 30), sigma=sigma,
                          alternative="greater", 
                          testType="paired", seed=1, pb=FALSE)

designObj2
#> 
#>  Savi Paired Sample T-Test Design
#> 
#>                       n1Plan, n2Plan = 30, 30
#> minimal standardised mean difference = 0.5656854
#>                          alternative = greater
#>                            power±2se = 0.719±0.0283044
#>                      e-variable type = mom
#>                      parameter: gMom = 0.16
#>                                alpha = 0.05
#>       Reject H_0: e-value >= 1/alpha = 20
#>              log(implied target)±2se = 3.882818±0.1586015
#> 
#> Timestamp: 2026-09-05 21:10:15 CEST

This reveals that due to budget constraints the experiment only has about 71.9% chance to detect the minimal clinical relevant effect size. If these chances are viewed as too slim, then we can either request more funds to invite more participants in the study, or prospectively decide that it is futile to conduct this experiment and spend our time and efforts on different endeavours instead.

Case (III): A desirable power and a budget constrained number of samples known. Target: Detectable effect size.

It is not always clear what the minimal clinically relevant effect size is, and provided with a budget for say, nMax=50, participants and a tolerable type I and desired power, we might be interested in finding the smallest detectable effect size. To do so, we run the following code:

# Recall:
# alpha <- 0.05
# power <- 0.8
designObj3 <- designSaviT(nPlan=c(50, 50), 
                          alpha=alpha, power=power,
                          sigma=sigma,
                          alternative="greater",
                          testType="paired")
designObj3
#> 
#>  Savi Paired Sample T-Test Design
#> 
#>                       n1Plan, n2Plan = 50, 50
#> minimal standardised mean difference = 0.5245391
#>                          alternative = greater
#>                      e-variable type = mom
#>                      parameter: gMom = 0.1375707
#>                                alpha = 0.05
#>       Reject H_0: e-value >= 1/alpha = 20
#> 
#> Timestamp: 2026-09-06 00:30:34 CEST
#> 
#> Note: The reported deltaMin is based on the batch analysis.

This shows that if we have budget for an experiment with 50 paired samples, and we analyse the data once after $n=50$, then we are able to detect a true effect size of 0.5245391 with 80% chance. This 0.5245391, however, is conservative, as when we monitor the e-value as the data come in, we will be able to detect even smaller effect sizes. If field experts believe that the true effect is smaller than 0.5245391, then we can again prospective decide against conducting this experiment, or we can obtain more funds.

2. Inference with savi tests: Optional stopping

In this section we highlight the point that the planned sample size for an e-value test is not rigid. There are roughly three cases: (1) the true effect size equals the minimal clinically relevant one, (2) the true effect size is larger than the minimal clinically relevant effect size, or (3) the true effect size is smaller than the minimal clinically relevant effect size. The main feature of e-variables is that we can monitor the e-value as the data come in, and stop the experiment early whenever the evidence is convincing. This occurs with high chance in scenarios (1) and (2). In scenario (3) the e-values at the planned sample size will typically be promising, but smaller than, say, $1/\alpha=20$. Continuing the experiment beyond the planned sample size does not hinder the validity of the e-value test, which we elaborate on in the next section on optional continuation. This and the next section thus shows that e-value tests are robust to both optional stopping and continuation, which implies that if the null hypothesis of no effect holds true, then there is less than $\alpha$ chance that the E-variable will ever reject the null.

We focus on scenarios (1) and (2) by first illustrating the operational characteristics of the savi test under the null, before we demonstrate its performance under the alternative.

Savi tests conserve the type I error rate: Batch analysis

We first show that the type I error is preserved for the batch analysis, that is, when the data are only analysed once at nPlan.

set.seed(1)
preData <- rnorm(n=designObj$nPlan[1], mean=120, sd=15)
postData <- rnorm(n=designObj$nPlan[2], mean=120, sd=15)
# Thus, deltaTrue=0
saviTTest(x=preData, y=postData, 
          designObj=designObj, paired=TRUE)
#> 
#>  Savi Paired Sample T-Test
#> 
#> data:  preData and postData. n1 = 36, n2 = 36
#> estimates: mean of the differences = -2.5799
#> 95 percent confidence sequence:
#>  -10.790209   7.631511
#> 
#> test: t = -0.78363, gMom = 0.16, type = mom
#> e-value = 0.020195 >= 1/alpha = 20 : FALSE
#> 
#> alternative hypothesis: true difference in means ('x' minus 'y') is greater than 0 
#> 
#> design: the test was designed with alpha = 0.05
#> for experiments with n1Plan = 36, n2Plan = 36
#> for minimal relevant standardised mean difference = 0.56569 (greater)

or equivalently with syntax closely resembling the standard t.test syntax:

savi.t.test(x=preData, y=postData, 
            designObj=designObj, paired=TRUE)
#> 
#>  Savi Paired Sample T-Test
#> 
#> data:  preData and postData. n1 = 36, n2 = 36
#> estimates: mean of the differences = -2.5799
#> 95 percent confidence sequence:
#>  -10.790209   7.631511
#> 
#> test: t = -0.78363, gMom = 0.16, type = mom
#> e-value = 0.020195 >= 1/alpha = 20 : FALSE
#> 
#> alternative hypothesis: true difference in means ('x' minus 'y') is greater than 0 
#> 
#> design: the test was designed with alpha = 0.05
#> for experiments with n1Plan = 36, n2Plan = 36
#> for minimal relevant standardised mean difference = 0.56569 (greater)

The following code replicates this simulation a 1,000 times and shows that in only a few cases will the E-variable cross the boundary of $1/\alpha$ under the null:

# alpha <- 0.05
nSim <- 1000

set.seed(1)
eValues <- replicate(n=nSim, expr={
  preData <- rnorm(n=designObj$nPlan[1], mean=120,
                   sd=15)
  postData <- rnorm(n=designObj$nPlan[2], mean=120,
                    sd=15)
  saviTTest(x=preData, y=postData, 
            designObj=designObj, paired=TRUE)$eValue}
)

mean(eValues >= 20)
#> [1] 0.003
mean(eValues >= 20) <= alpha
#> [1] TRUE

Hence, in this simulation with the null hypothesis holding true and if the savi test is only conducted once at the planned sample size, then in 3 out of 1000 experiments the null hypothesis was falsely rejected.

Savi tests allow for early stopping without inflating the type I error rate above the tolerable $\alpha$-level

What makes the savi tests in this package particularly interesting is that they allow for early stopping without the test ever exceeding the tolerable type I error rate of $\alpha$. This means that the e-value can be monitored as the data come in, and when there is a sufficient amount of evidence against the null, that is, whenever the e-value is larger than $ 1/$, the experiment can be stopped early. This puts fewer patients at risk, and allows for more efficient scientific scrutiny.

Note that not all E-variables necessarily allow for optional stopping: this only holds for some special E-variables, that are also test martingales. More information can be found, for example, in Rosanne’s master thesis, Chapter 5.

Optional stopping is problematic for p-values

Optionally stopping will never cause the savi test to over-reject the null hypothesis, whereas tracking the classical p-value tests and rejecting the null whenever it dips below, say, $\alpha=0.05$ will lead to an over-inflation of the type I error. In other words, optional stopping with these p-value tests leads to an increased risk of falsely claiming that a medicine is effective, while in reality it is not.

The following code replicates 1000 experiments and each data set is generated with a true effect size set to zero. We first collect the t-statistics across the 1000 data sets, and time/sample size 36

nSim <- 1000
muGlobal <- 120
n1 <- designObj$nPlan[1]

nullData <- generateNormalData(
  designObj$nPlan, muGlobal=muGlobal,
  nSim=nSim, deltaTrue=0, seed=1,
  sigma=sigma)

# Used to vectorise the computations of the t-statistic
n1Vector <- 1:n1
nuVector <- n1Vector-1

# Here we store all the t statistics across the 
# number of simulations (nSim) and time (n1)
tMatrix <- matrix(nrow=nSim, ncol=n1)

for (sim in 1:nSim) {
  dataGroup1 <- nullData$dataGroup1[sim, ]
  dataGroup2 <- nullData$dataGroup2[sim, ]
  
  differenceScore <- dataGroup1-dataGroup2
  
  # Vector of mean differences
  meanDiffVector <- 
    1/n1Vector*cumsum(differenceScore)
  
  # Vector of standard deviations
  sdMeanDiff <- sqrt(
    1/nuVector*(cumsum(differenceScore^2)-n1Vector*meanDiffVector^2)
  )
  
  tMatrix[sim, ] <- 
    sqrt(n1Vector)*meanDiffVector/sdMeanDiff
}

# The first t-value is undefined
tMatrix[, 1] <- 0

For each data set a sequential analysis is run. We will see that we quickly exceed the tolerable type I error if we monitor the p-value as the data come in, and reject the null as soon as the p-value dips below $\alpha=0.05$.

# Here we store all the p-values across the 
# number of simulations (nSim) and time (n1)
allPValues <- matrix(nrow=nSim, ncol=n1)

# Whenever this vector has a 1 it indicates that the simulate data 
# yielded a "significant" p-value, despite the data being generated under the null
pValueUnderAlpha <- vector("integer", length=nSim)

# This indicates the first time an experiment yielded a p-value < alpha=0.05
# Default is Inf, which indicates that the p-value dip below alpha
firstPassageTime <- rep(Inf, times=nSim)

for (sim in 1:nSim) {
  tVector <- tMatrix[sim, ]
  
  for (i in 2:n1) {
    currentPValue <- 1-stats::pt(tVector[i], df=nuVector[i])
    allPValues[sim, i] <- currentPValue
    
    if (currentPValue < alpha && 
        pValueUnderAlpha[sim]!=1) {
      pValueUnderAlpha[sim] <- 1
      firstPassageTime[sim] <- i
      break()
    }
  }
}
  
numberOfDippingExperimentsAtTimeN <- integer(n1)

for (i in 1:n1) {
  numberOfDippingExperimentsAtTimeN[i] <- 
    sum(firstPassageTime <= i)
}

pValueFalseRejects <- numberOfDippingExperimentsAtTimeN/nSim

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
plot(1:n1, 100*pValueFalseRejects, type="l", xlab="n",
     ylab="Type I error (%)", ylim=c(0, 25),
     lwd=2, col=freqColours[1])
lines(c(1, n1), c(5, 5), lwd=2, lty=2)

After the second evaluation a total of 93 out of nSim=1000 experiments yielded a significant p-value after the first or second observation, that is, 9.3, which is already larger than the tolerable 5%. This emphasises the point that the significant test $p < 0.05$ should be conducted once, and only once.

Optional stopping does not cause savi tests to over-reject the null

We repeat the simulation, but now with e-values instead. With e-values we can reject the null and stop the experiment as soon as the e-value exceeds $1 / \alpha = 20$. Importantly, this procedure will not yield more than the tolerable 5% false positive errors up to nPlan, but also beyond, that is, if we add more participants, see the sections concerned with optional continuation.

eValueFalseRejects

# Here we store all the e-values across the 
# number of simulations (nSim) and time (n1)
eValues <- matrix(nrow=nSim, ncol=n1)

# This indicates whether a simulation yielded e >= 1/alpha
eOver <- vector("integer", length=nSim)

# This indicates the first time an experiment yielded e >= 1/alpha
# Default is Inf, which indicates that the e didn't cross 1/alpha
firstPassageTimeE <- rep(Inf, times=nSim)

# This is the e-value at the end time, or whenever 
# it exceeds the threshold of 1/alpha
eStopped <- numeric(nSim)

for (sim in 1:nSim) {
  tVector <- tMatrix[sim, ]
  
  for (i in 1:n1) {
    currentEValue <- saviTTestStat(
      tVector[i], parameter=designObj$parameter, 
      n1=n1Vector[i], n2=n1Vector[i], 
      paired=TRUE,  eType=designObj$eType)$eValue
    
    eValues[sim, i] <- currentEValue
    
    if (currentEValue >= 1/alpha && eOver[sim]!=1) {
      eOver[sim] <- 1
      firstPassageTimeE[sim] <- i
      eStopped[sim] <- currentEValue
    }
    
    if (i==n1 && eOver[sim]!=1) {
      eStopped[sim] <- currentEValue
    }
  }
}

trackCrossing <- integer(n1)

for (i in 1:n1) {
  trackCrossing[i] <- 
    sum(firstPassageTimeE <= i)
}

eValueFalseRejects <- trackCrossing/nSim

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
plot(1:n1, 100*eValueFalseRejects, type="l", 
     xlab="n", ylab="Type I error (%)", lwd=2,
     col=eColours[1], ylim=c(0, 5))
lines(c(1, n1), c(5, 5), lwd=2, lty=2)

Note that optional stopping always increases the chance of observing a false detection. For the savi test this increased to 2.8%, which is still below the tolerable 5%. On the other hand, tracking the p-value and rejecting the null as soon it falls below $\alpha$ leads to 25.4%, which is well above 5%.

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
plot(1:n1, 100*pValueFalseRejects, type="l", xlab="n",
     ylab="Type I error (%)",
     lwd=2, col=freqColours[1], ylim=c(0, 25))
lines(1:n1, 100*eValueFalseRejects, lwd=2,
     col=eColours[1])
lines(c(1, n1), c(5, 5), lwd=2, lty=2)

Savi tests detect the effect early if it is present: deltaTrue equal to deltaMin

In this section we illustrate the operational characteristics of the savi t-test under optional stopping, when the effect is present. The following code replicates 1000 experiments and each data set is generated with a true mean difference that equals the minimal clinically relevant effect size of $\delta_{\min}=0.57$. If the e-value does not exceed $1 / \alpha$, the experiment is run until all samples are collected as planned. The simulation described above can be run with the following command:

simDeltaTrueIsDeltaMin <- 
  sampleStoppingTimesSaviT(
    deltaTrue=deltaMin, alternative="greater", 
    testType="paired", sigma=sigma, 
    nMax=designObj$nPlan, seed=1,
    parameter=designObj$parameter,nSim=nSim)

mean(
  simDeltaTrueIsDeltaMin$eValuesStopped >=
    20)
#> [1] 0.801

The simulations confirm that there is an 80% chance of detecting the minimal clinically relevant effect if we run the test until the planned sample size. Typically, there is a discrepancy due to sampling error, which vanishes as the number of simulations increases.

To compute the mean and plot the distributions of the stopping times, we run the following code:

stoppingTimes <- simDeltaTrueIsDeltaMin$stoppingTimes

mean(stoppingTimes)

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes()
hist(stoppingTimes, 
     breaks=min(stoppingTimes):max(stoppingTimes), 
     xlim=c(0, designObj$nPlanBatch[1]),col=eColours[2],
     border=eColours[1], lwd=2, main="")

The histogram shows the full distribution of the times at which the experiment is stopped. For instance, 412 out of the 1000 experiments stopped before half the planned sample size. In these cases we were lucky and the effect was detected early. The last bar collects all experiments that ran until the planned sample sizes, thus, also those that did not lead to a null rejection at n=36. To see the distributions of stopping times of only the experiments where the null is rejected, we run the following code:

firstPassageTimeAltEqual <- stoppingTimes
firstPassageTimeAltEqual[
  which(as.integer(simDeltaTrueIsDeltaMin$breakVector)==1)] <- Inf

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes()
hist(firstPassageTimeAltEqual, 
     breaks=min(firstPassageTimeAltEqual):n1, 
     xlim=c(0, n1),col=eColours[2],
     border=eColours[1], lwd=2, main="")

This can also be visualised as follows:

trackCrossingAltEqual <- integer(n1)

for (i in 1:n1) {
  trackCrossingAltEqual[i] <- sum(firstPassageTimeAltEqual <= i)
}

eReject <- trackCrossingAltEqual/nSim

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
plot(1:n1, 100*eReject, type="l", 
     xlab="n", ylab="Correct rejections (%)", lwd=2,
     col=eColoursAlt[1], ylim=c(0, 80))
lines(c(1, n1), c(80, 80), lwd=2, lty=2)

Here the horizontal line represents the 80% targeted power, which is indeed reached at n=36 under optional stopping.

Savi tests detect the effect even earlier if deltaTrue is larger than deltaMin

What we believe is clinically minimally relevant might not match reality. One advantage of savi tests is that they perform even better, if the true effect size is larger than the minimal clinical effect size. This is illustrated with the following code

simDeltaTrueLargerDeltaMin <- 
  sampleStoppingTimesSaviT(
    deltaTrue=1.2*deltaMin, alternative="greater",
    testType="paired", sigma=sigma, 
    nMax=designObj$nPlan, seed=1,
    parameter=designObj$parameter,nSim=nSim)

mean(
  simDeltaTrueLargerDeltaMin$eValuesStopped >=
    20)
#> [1] 0.935

With a larger true effect size, the power increased to 93.5%. More importantly, this increase is picked up earlier by the designed savi test, and optional stopping allows us to act on this. Note that the average stopping time is now decreased, from 21.672 to 17.517. This is apparent from the fact that the histogram of stopping times is now shifted to the left:

stoppingTimes <- 
  simDeltaTrueLargerDeltaMin$stoppingTimes

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes()
hist(stoppingTimes, 
     breaks=min(stoppingTimes):max(stoppingTimes), 
     xlim=c(0, designObj$nPlanBatch[1]),col=eColours[2],
     border=eColours[1], lwd=2, main="")

Hence, this means that if the true effect is larger than what was planned for, the savi test will detect this larger effect earlier on, which results in a further increase of efficiency, as shown by the following plot.

firstPassageTimeAltLarger <- 
  simDeltaTrueLargerDeltaMin$stoppingTimes
firstPassageTimeAltLarger[
  which(as.integer(simDeltaTrueLargerDeltaMin$breakVector)==1)] <- Inf

trackCrossingAltLarger <- integer(n1)

for (i in 1:n1) {
  trackCrossingAltLarger[i] <- sum(firstPassageTimeAltLarger <= i)
}

eReject <- trackCrossingAltLarger/nSim

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
plot(1:n1, 100*eReject, type="l", 
     xlab="n", ylab="Correct rejections (%)", lwd=2,
     col=eColoursAlt[1], ylim=c(0, 90))
lines(c(1, n1), c(80, 80), lwd=2, lty=2)

The plot visualises the increase in power for a true effect size larger than the minimal clinically relevant one.

The last scenario with deltaTrue smaller than the minimally clinically relevant effect size is discussed in the context of optional continuation.

3. Optional Continuation

Savi tests, as we will show below, conserve the type I error rate under optional continuation. Optional continuation implies gathering more samples than was planned for because, for instance, (1) more funding became available and the experimenter wants to learn more, (2) the evidence looked promising, (3) a reviewer or editor urged the experimenter to collect more data, or (4) other researchers attempt to replicate the first finding.

A natural way to deal with the first three cases is by computing an e-value over the combined data set. This is permitted if the data come from the same population, and if the E-variable used is a test martingale, which is the case for the problem at hand.

Replication attempts, however, are typically based on samples from a different population. One way to deal with this is by multiplying the e-value computed from the original study with the e-value computed from the replication attempt. In this situation, the e-value formula for the replication study could also be redesigned through the function, for example when more information on nuisance parameters or effect size has become available to design a more powerful test.

We show that both procedures are safe, that is, neither results in the e-value test over-rejecting the null, whenever it holds true, as is the case with classical p-values. We first show that optional continuation with p-values is problematic.

Optional continuation is problematic for p-values and leads to overinflating the type I error rate

Firstly, we show that optional continuation also causes p-values to over-reject the null. In the previous section we saw that optional stopping causes the p-value to falsely reject the null with about 20% chance, much higher than the tolerable 5%. We consider the situation where the p-value is performed once, at the end of the trial, at n1. The non-significant studies then get extended with a second batch of data. We will see that this procedure of selectively continuing non-significant experiments causes the collective rate of false null rejections to be larger than $\alpha$.

We begin with correctly computed p-value analyses at the last sample n1 based on the following t-statistics.

tAtN1 <- numeric(nSim)

for (sim in 1:nSim) {
  dataGroup1 <- nullData$dataGroup1[sim, ]
  dataGroup2 <- nullData$dataGroup2[sim, ]
  
  meanDiff <- 
    mean(dataGroup1-dataGroup2)
  
  sdMeanDiff <- sd(dataGroup1-dataGroup2)
  
  tAtN1[sim] <- sqrt(n1)*meanDiff/sdMeanDiff
}

The p-values after the first batch of data are computed as follows:

pValuesBatch1 <- numeric(nSim)

for (i in 1:nSim) {
  pValuesBatch1[i] <- 1-stats::pt(tAtN1[i], df=n1-1)
}

mean(pValuesBatch1 < alpha)
#> [1] 0.046

Hence, after a first batch of data, we get 46 incorrect null rejections out of 1000 experiments (4.6%).

The following code continues only the non-significant 954 experiments with a second batch of data, all also generated under the null.

nullData2 <- generateNormalData(
  designObj$nPlan, muGlobal=muGlobal,
  nSim=nSim, deltaTrue=0, seed=2,
  sigma=sigma)

tAtN2 <- numeric(nSim)

for (sim in 1:nSim) {
  dataGroup1 <- c(nullData$dataGroup1[sim, ], 
                  nullData2$dataGroup1[sim, ])
  dataGroup2 <- c(nullData$dataGroup2[sim, ], 
                  nullData2$dataGroup2[sim, ])
  
  meanDiff <- mean(dataGroup1-dataGroup2)
  sdMeanDiff <- sd(dataGroup1-dataGroup2)
  
  tAtN2[sim] <- sqrt(2*n1)*meanDiff/sdMeanDiff
}

rejectedIndeces <- which(pValuesBatch1 < alpha)
notRejectedIndeces <- which(pValuesBatch1 >= alpha)

pValuesBatch2 <- numeric(nSim)
pValuesBatch2[rejectedIndeces] <- 
  pValuesBatch1[rejectedIndeces]

for (j in notRejectedIndeces) {
  pValuesBatch2[j] <- 1-stats::pt(tAtN2[j], df=2*n1-1)
}

mean(pValuesBatch2[notRejectedIndeces] < alpha)
#> [1] 0.02830189

By selectively extending the non-significant results of the first batch with a second batch of data, we got an additional 27 false rejection. This brings the collective total to 73 of false rejections out of a total of 1000 studies. That is, a false positive rate of 7.3%, which is above the tolerable 5%.

The reason why p-values over-reject the null under optional stopping and optional continuation is due to p-values being uniformly distributed under the null.

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
hist(pValuesBatch1, col=freqColours[1])
abline(v=0.05)

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
hist(pValuesBatch2[notRejectedIndeces], col=freqColours[2])
abline(v=0.05)

As such, if the null holds true and the number of samples increases, then the p-value meanders between 0 and 1, thus, eventually crossing any fixed $\alpha$-level. This simulation reiterates the point that classical p-value tests only retain their type I error control if they are used once, and only once.

i. Optional continuation by extending the experiment does not result in savi tests exceeding the tolerable $\alpha$-level

In this section we provide more insight to what it means for e-variables to be anytime-valid. More precisely, we show that type I error control is retained even if we keep on monitoring the e-value after the planned sample size. We do so in a more difficult setting compared to the previous p-value analysis. Instead of analysing the data in batches as we did with the p-value, let us illustrate this by extending the studies with continuous streams of data, say ten times the planned sample sizes. That is, we show that type I error is retained even after 360 looks, and recall that the type I error guarantee of the p-value breaks after two looks within a batch, or after two batches. The following code is used to illustrate this sample size independent type I error control.

# Additional data under the null
nullData3 <- generateNormalData(
  9*designObj$nPlan, muGlobal=muGlobal,
  nSim=nSim, deltaTrue=0, seed=3,
  sigma=sigma)

# All the z statistics across the 
# number of simulations (nSim) and time (10*n1)
tMatrixAll <- matrix(nrow=nSim, ncol=10*n1)

# Used to vectorise the computations for the the z-statistic
n1Vector <- 1:(10*n1)
nuVector <- n1Vector-1

for (sim in 1:nSim) {
  dataGroup1 <- c(nullData$dataGroup1[sim, ],
                  nullData3$dataGroup1[sim, ])
  dataGroup2 <- c(nullData$dataGroup2[sim, ],
                  nullData3$dataGroup2[sim, ])
  
  differenceScore <- dataGroup1-dataGroup2
  
  meanDiffVector <- 
    1/n1Vector*cumsum(differenceScore)
  
  # Vector of standard deviations
  sdMeanDiff <- sqrt(
    1/nuVector*(cumsum(differenceScore^2)-n1Vector*meanDiffVector^2)
  )
  
  tMatrixAll[sim, ] <- 
    sqrt(n1Vector)*meanDiffVector/sdMeanDiff
}

tMatrixAll[, 1] <- 0

# Here we store all the e-values across the 
# number of simulations (nSim) and time (n1)
allEValues <- matrix(nrow=nSim, ncol=10*n1)
allEValues[, 1:n1] <- eValues

eOverOptioCont <- eOver

for (sim in 1:nSim) {
  tVector <- tMatrixAll[sim, ]
  
  for (i in (n1+1):(10*n1)) {
    currentEValue <- saviTTestStat(
      tVector[i], parameter=designObj$parameter, 
      n1=n1Vector[i], n2=n1Vector[i], 
      paired=TRUE,  sigma=sigma,
      eType=designObj$eType)$eValue
    
    allEValues[sim, i] <- currentEValue
    
    if (currentEValue >= 1/alpha && eOverOptioCont[sim]!=1) {
      eOverOptioCont[sim] <- 1
      firstPassageTimeE[sim] <- i
    }
  }
}

trackCrossingOptioCont <- integer(10*n1)

for (i in 1:(10*n1)) {
  trackCrossingOptioCont[i] <- 
    sum(firstPassageTimeE <= i)
}

eValueFalseRejects2 <- 
  trackCrossingOptioCont/nSim

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
plot(1:(10*n1), 100*eValueFalseRejects2, type="l", 
     xlab="n", ylab="Type I error (%)", lwd=2,
     col=eColours[1], ylim=c(0, 5))
lines(c(1, 10*n1), c(5, 5), lwd=2, lty=2)

The original paper on e-values by Grunwald, de Heide and Koolen and Howard, Ramdas, McAuliffe, and Sekhon provide mathematical proofs showing that the type I error will never exceed the tolerable type I error rate.

The simulations show that the sample size computed by the design function is indeed a non-rigid planned sample size. The planned sample size is actually only concerned with the alternative and does not involve the null. This is in contrast to certain alpha-spending procedures that only provide type I error up to a maximum sample size, as all alpha is spent at that point.

Under the null, the behaviour of e-variables is very different to that of p-values, which meander between zero and one. In contrast, e-variables slowly drift to zero under the null. This is equivalent to drifting towards -infinity on the logarithmic scale. The following plot illustrates this slow drift to zero.

lowerQuartileLogEValueNull <- numeric(10*n1)
medianLogEValueNull <- numeric(10*n1)
upperQuartileLogEValueNull <- numeric(10*n1)

for (j in 1:(10*n1)) {
  brie <- quantile(log(allEValues[, j]))
  lowerQuartileLogEValueNull[j] <- brie[2]
  medianLogEValueNull[j] <- brie[3]
  upperQuartileLogEValueNull[j] <- brie[4]
}

nDomain <- 1:(10*n1)

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
plot(nDomain, medianLogEValueNull, col="black", lwd=2, 
     ylim=c(-6, 3), type="l", xlab="n", 
     ylab="log(eValues)")
lines(nDomain, lowerQuartileLogEValueNull, col=eColours[1],
      lwd=2, lty=1)
lines(nDomain, upperQuartileLogEValueNull, col=eColours[1],
      lwd=2, lty=1)
lines(c(0, 10*n1), c(log(20), log(20)), lwd=2, col="grey", lty=2)

The black curve depicts the median of the e-variable and the blue curves represent the 25% and 75% percentile of the e-variable distribution under the null. The horizontal grey line depicts $\log(1/\alpha) \approx 3$ for $\alpha=0.05$. The plot shows that under the null it gets increasingly hard for e-variables to cross the threshold of $\log(1/\alpha)$ as the sample size increase. This also illustrates why the e-variable’s marginal increase in type I error slowly diminishes.

When the effect is present optional continuation results in savi tests correctly rejecting the null

The slow drift of the sampling distribution of e-values to smaller values is replaced by a fast drift to large values whenever there is an effect. We consider the situation where the study is continued after the planned sample size, but with deltaTrue = 0.5, thus, smaller than deltaMin = 0.57.

# Data under the alternative
altData <- generateNormalData(
  10*designObj$nPlan, muGlobal=muGlobal,
  nSim=nSim, deltaTrue=0.5, seed=2,
  sigma=sigma)

# All the z statistics across the 
# number of simulations (nSim) and time (10*n1)
tMatrixAllAlt <- matrix(nrow=nSim, ncol=10*n1)

# Used to vectorise the computations for the z-statistic
n1Vector <- 1:(10*n1)

for (sim in 1:nSim) {
  dataGroup1 <- altData$dataGroup1[sim, ]
  dataGroup2 <- altData$dataGroup2[sim, ]
  
  meanDiffVector <- 
    1/n1Vector*cumsum(dataGroup1-dataGroup2)
  
  # The variance of the sum x + (-y) is the sum of the two variances
  # Thus, 2*sigma^2
  sdMeanDiff <- sqrt(2)*sigma
  
  tMatrixAllAlt[sim, ] <- 
    sqrt(n1Vector)*meanDiffVector/sdMeanDiff
}


# Here we store all the e-values across the 
# number of simulations (nSim) and time (n1)
allEValuesAlt <- matrix(nrow=nSim, ncol=10*n1)

eOverAlt <- integer(nSim)
firstPassageTimeEAlt <- rep(Inf, nSim)
eStoppedAlt <- numeric(nSim)

for (sim in 1:nSim) {
  tVector <- tMatrixAllAlt[sim, ]
  
  for (i in 1:(10*n1)) {
    currentEValue <- saviZTestStat(
      tVector[i], parameter=designObj$parameter, 
      n1=n1Vector[i], n2=n1Vector[i], 
      paired=TRUE,  sigma=sigma,
      eType=designObj$eType)$eValue
    
    allEValuesAlt[sim, i] <- currentEValue
    
    if (currentEValue >= 1/alpha && eOverAlt[sim]!=1) {
      eOverAlt[sim] <- 1
      firstPassageTimeEAlt[sim] <- i
      eStoppedAlt[sim] <- currentEValue
    }
    
    if (i==n1 && eOverAlt[sim]!=1) {
      eStoppedAlt[sim] <- currentEValue
    }
  }
}

trackCrossingOptioCont <- integer(10*n1)

for (i in 1:(10*n1)) {
  trackCrossingOptioCont[i] <- 
    sum(firstPassageTimeEAlt <= i)
}

eValueCorrectRejects <- 
  trackCrossingOptioCont/nSim

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
plot(1:(10*n1), 100*eValueCorrectRejects, type="l", 
     xlab="n", ylab="Correct rejections (%)", lwd=2,
     col=eColoursAlt[1])#, #ylim=c(0, 5))
lines(c(1, 10*n1), c(5, 5), lwd=2, lty=2)

The plot illustrates that for the smaller deltaTrue, there is only about 30.6% power to correctly reject the null at the planned sample size of nPlan=36. Since it is not possible to over-reject the null with an e-variable under optional continuation, we can simply continue the study whenever the evidence seems convincing. Doing so thus allows effects to be detected even if they are smaller than anticipated.

The mentioned fast drift towards higher values is illustrated by the following plot:

lowerQuartileLogEValueAlt <- numeric(10*n1)
medianLogEValueAlt <- numeric(10*n1)
upperQuartileLogEValueAlt <- numeric(10*n1)

for (j in 1:(10*n1)) {
  brie <- quantile(log(allEValuesAlt[, j]))
  lowerQuartileLogEValueAlt[j] <- brie[2]
  medianLogEValueAlt[j] <- brie[3]
  upperQuartileLogEValueAlt[j] <- brie[4]
}

nDomain <- 1:(10*n1)

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
plot(nDomain, medianLogEValueAlt, col="red", lwd=2, 
     ylim=c(-6, 15), type="l", xlab="n", 
     ylab="log(eValues)")
lines(c(0, 10*n1), c(log(20), log(20)), lwd=2, col="grey", lty=2)
lines(nDomain, lowerQuartileLogEValueAlt, col=eColoursAlt[1],
      lwd=2, lty=1)
lines(nDomain, upperQuartileLogEValueAlt, col=eColoursAlt[1],
      lwd=2, lty=1)

lines(nDomain, medianLogEValueNull, col="black", lwd=2)
lines(nDomain, lowerQuartileLogEValueNull, col=eColours[1],
      lwd=2, lty=1)
lines(nDomain, upperQuartileLogEValueNull, col=eColours[1],
      lwd=2, lty=1)

The red curve depicts the median of the e-variable and the brown curves represent the 25% and 75% percentile of the e-variable distribution under the alternative. The horizontal grey line depicts $\log(1/\alpha) \approx 3$ for $\alpha=0.05$. The plot shows that under the alternative it gets increasingly easier for e-variables to cross the threshold of $\log(1/\alpha)$ as the sample size increases.

4. Optional contination and meta-analysis: Extending the research endeavours through replication studies

It is not always appropriate to compute e-values over combined data sets, in particular for replication attempts where the original experiment is performed on a different population. Instead of combining the data we combine the evidence of the individual studies by multiplying their respective e-values. This procedure is also safe under optional continuation, as the type I error is also conserved.

Type I error is guaranteed when multiplying e-values

To demonstrate type I error control in a meta-analytical setting we consider the data generated in the optional stopping section as original studies, and follow-up studies also with no effect. With the aforementioned original studies we carry over a false positive rate of 2.8%. As with the original study we will also continuously monitor the e-values in the follow-up studies. By multiplying e-values we might gain efficiency as follows: If an original study resulted in an e-value of 10 at nPlan, then we only need an e-value of 2 in the follow-up study as the data come in for the product to be higher than the threshold of 20 when $\alpha=0.05$.

As in the original study the data are generated under the null, but we assume that the drug is now administered to a clinical group that has a lower overall baseline blood pressure of $\mu_{g}=90$ mmHg and standard deviation of $\sigma=6$. We assume that the follow-up study has a planned sample size of 2 times nPlan, but the statement also holds true for follow-up studies with a smaller sample size. The following code selectively continues the original studies and save them in the list rep2.

rep2 <- selectivelyContinueZOrTTestData(
  designObj, n1New=2*n1, testName="T-Test",
  muGlobal=90, sigma=6, deltaTrue=0, nSim=nSim,
  eValuesOld=eValues, eOverOld=eOver, 
  trackCrossingOld=trackCrossing, 
  firstPassageTimeOld=firstPassageTimeE,
  eStoppedOld=eStopped, seed=2)

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
plot(1:(3*n1), 100*rep2$trackCrossing/nSim, type="l", 
     xlab="n", ylab="Type I error (%)", lwd=2,
     col=eColours[1], ylim=c(0, 5))
lines(c(1, 3*n1), c(5, 5), lwd=2, lty=2)

rep2$extraRejections
#> [1] 8

By selectively continuing the not rejected studies we now incurred an additional 8 false discoveries (3.6% error).

Let’s consider another selective replication round with a planned sample size about half nPlan, but with the drug administered to yet another population.

rep3 <- selectivelyContinueZOrTTestData(
  designObj, n1New=ceiling(0.48*n1), testName="T-Test",
  muGlobal=100, sigma=8, deltaTrue=0, nSim=1000,
  eValuesOld=rep2$eValues, eOverOld=rep2$eOver, 
  trackCrossingOld=rep2$trackCrossing, 
  firstPassageTimeOld=rep2$firstPassageTime,
  eStoppedOld=rep2$eStopped, seed=3)

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
plot(1:(length(rep3$trackCrossing)), 100*rep3$trackCrossing/nSim, type="l", 
     xlab="n", ylab="Type I error (%)", lwd=2,
     col=eColours[1], ylim=c(0, 5))
lines(c(1, length(rep3$trackCrossing)), c(5, 5), lwd=2, lty=2)

rep3$extraRejections
#> [1] 0

By selectively continuing the not rejected studies we now incurred no additional 0 false discoveries retaining the cumulative 3.6% type I error rate.

A fourth replication round with

rep4 <- selectivelyContinueZOrTTestData(
  designObj, n1New=ceiling(3.5*n1), testName="T-Test",
  muGlobal=150, sigma=19, deltaTrue=0, nSim=1000,
  eValuesOld=rep3$eValues, eOverOld=rep3$eOver, 
  trackCrossingOld=rep3$trackCrossing, 
  firstPassageTimeOld=rep3$firstPassageTime,
  eStoppedOld=rep3$eStopped, seed=4)

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
plot(1:(length(rep4$trackCrossing)), 100*rep4$trackCrossing/nSim, type="l", 
     xlab="n", ylab="Type I error (%)", lwd=2,
     col=eColours[1], ylim=c(0, 5))
lines(c(1, length(rep4$trackCrossing)), c(5, 5), lwd=2, lty=2)

rep4$extraRejections
#> [1] 2

Again we did not incur additional false discoveries, and the process can be repeated ad nauseam and the type I error will always remain under the tolerable $\alpha=0.05$. The reason for this is again the slowly drift of e-variables towards 0 under the null, thus, -infinity on the logarithmic scale:

totalN <- length(rep4$trackCrossing)

lowerQuartileLogEValueMetaNull <- numeric(totalN)
medianLogEValueMetaNull <- numeric(totalN)
upperQuartileLogEValueMetaNull <- numeric(totalN)

for (i in 1:totalN)  {
  brie <- quantile(log(rep4$eValues[, i]))
  
  lowerQuartileLogEValueMetaNull[i] <- brie[2]
  medianLogEValueMetaNull[i] <- brie[3]
  upperQuartileLogEValueMetaNull[i] <- brie[4]
}

yMin <- floor(min(lowerQuartileLogEValueMetaNull))

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
plot(1:totalN, medianLogEValueMetaNull, col="black", lwd=2, 
     ylim=c(yMin, 3), type="l", xlab="n", 
     ylab="log(eValues)")
lines(1:totalN, lowerQuartileLogEValueMetaNull, col=eColours[1],
      lwd=2, lty=1)
lines(1:totalN, upperQuartileLogEValueMetaNull, col=eColours[1],
      lwd=2, lty=1)
lines(c(0, totalN), c(log(20), log(20)), lwd=2, col="grey", lty=2)
lines(c(n1, n1), c(yMin, 3), lty=2, col="lightgrey")
lines(c(3*n1, 3*n1), c(yMin, 3), lty=2, col="lightgrey")
lines(c(3.5*n1, 3.5*n1), c(yMin, 3), lty=2, col="lightgrey")

The vertical lines indicate where a follow-up study was initiated.

Meta-analysis e-values under the alternative

As original experiments we now take the e-values from the optional stopping simulation study with deltaTrue equal to deltaMin. The code below selectively continues studies where the null was not rejected. We replicate the study in a population with different nuisance parameters, e.g. $\mu_{g}=110$ and $\sigma=50$, thus, much more spread out than in the original studies. The following code illustrates the correct rejection rate when the effect size of the replication is similar to that of the original study.

repAlt <- selectivelyContinueZOrTTestData(
  designObj, n1New=ceiling(2*n1), testName="T-Test",
  muGlobal=145, sigma=15, deltaTrue=deltaMin, nSim=1000,
  eValuesOld=simDeltaTrueIsDeltaMin$samplePaths, 
  eOverOld=!simDeltaTrueIsDeltaMin$breakVector, 
  trackCrossingOld=trackCrossingAltEqual, 
  firstPassageTimeOld=firstPassageTimeAltEqual,
  eStoppedOld=simDeltaTrueIsDeltaMin$eValuesStopped, seed=6)

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
plot(1:(length(repAlt$trackCrossing)), 100*repAlt$trackCrossing/nSim, type="l", 
     xlab="n", ylab="Correct rejections (%)", lwd=2,
     col=eColoursAlt[1], ylim=c(0, 100))
lines(c(1, length(repAlt$trackCrossing)), c(80, 80), lwd=2, lty=2)
lines(c(n1, n1), c(0, 100), lty=2, col="lightgrey")

The vertical grey line is drawn at the nPlan of the original study at which we correctly reject the null with 80.1% power. On the right of this grey line we see that the number of correct rejections further increases by combining the original e-values with e-values from replication attempts with different nuisance parameters.

Conclusion

We believe that optional continuation is essential for (scientific) learning, as it allows us to revisit uncertain decisions such as ($p < \alpha$ and e-value $\geq 1/\alpha$) either by extending an experiment directly, or via replication studies in a meta-analysis. Hence, we view learning as an ongoing process, which requires that inference becomes more precise as data accumulate. The inability of p-values to conserve the $\alpha$-level under optional continuation, however, is at odds with this view –by gathering more data after an initial look, the inference becomes less precise, as the likelihood of the null being true after observing $p < \alpha$ increases beyond what is tolerable.

Savi tests on the other hand benefit from more data, as the chance of seeing e-value $\geq 1/\alpha$ (slowly) decreases when the null is true, whereas it (quickly) increases when the alternative is true, as the number of samples increases.

Flexible Safe Anytime-Valid Tests for Practical Scenarios

Alexander Ly, Udo Boehm and Rosanne Turner