---
title: "Flexible Safe Anytime-Valid Tests for Practical Scenarios"
author: "Alexander Ly, Udo Boehm and Rosanne Turner"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Flexible Safe Anytime-Valid Tests for Practical Scenarios}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.height = 4,
  fig.width = 8
)
```

_Safe anytime-valid inference_ (savi) is a collective name for a new method of inference based on _e-values_ (instead of _p-values_). The original paper on e-values by Grunwald, de Heide and Koolen can be found [on the website of Journal of the Royal Statistical Society Series B: Statistical Methodology](https://doi.org/10.1093/jrsssb/qkae069). Another good source on the mathematical theory behind savi tests and e-values is provided by [Howard, Ramdas, McAuliffe, and Sekhon](https://doi.org/10.1214/20-AOS1991), the book ["Hypothesis testing with e-values"](https://doi.org/10.1561/3600000002) by Ramdas and Wang, also available on [arXiv:2410.23614](https://arxiv.org/abs/2410.23614), and the tutorial paper 
["A Tutorial on Safe Anytime-Valid Inference: Practical Maximally Flexible Sampling Designs for Experiments Based on e-Values"](https://osf.io/preprints/psyarxiv/h5vae) by Ly, Boehm, Grunwald, Ramdas, and Van Ravenzwaaij. 

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. 

>Technically, E-variables, are non-negative random variables (test statistics) that have an expected value of at most one under the null hypothesis. The E-variable can be interpreted as an gamble against the null hypothesis in which an investment of 1\$ returns E\$ whenever the null hypothesis fails to hold true. Hence, the larger the observed e-value, the larger the incentive to reject the null (see the [original paper](https://arxiv.org/abs/1906.07801)). 

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](#designT). 

Secondly, simulations are run to show that savi tests indeed conserve the type I error guarantee under optional stopping [for testing means](#optStopT). 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](#optContT).

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. <!--, and flexibility is required to act on new information.Furthermore, the failure of p-values to be robust to optional continuation implies that one cannot build on previously acquired data, thus, engage in (scientific) learning. 

simulations are run to show that savi tests also conserve the type I error guarantee under optional continuation. This implies that 

we show that the behaviour of savi tests under optional continuation. 
-->

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`:

```{r install, eval=FALSE}
install.packages("safestats")
```

The development version can be found on [GitHub](https://github.com/alexanderlynl/safestats), which can be installed with the `remotes` package from [CRAN](https://cran.r-project.org/package=devtools) by entering in `R`:

```{r remotes, eval=FALSE}
remotes::install_github("AlexanderLyNL/safestats", build_vignettes = TRUE)
```

The command 

```{r setup}
library(safestats)
```

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

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

# 1. <a name = "designT">Designing savi t-test experiments</a>
## 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. <!--This can only be accomplished by tests that never reject the null. -->
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.  <!-- In other words, whenever there is an effect, then 80 out of a 100 experiment should reject the null hypothesis of no effect. -->

## 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. <!--The $\sqrt{2}$-term in the denominator is a result of the measurements being paired. -->

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: 

```{r}
alpha <- 0.05
power <- 0.8
deltaMin <- 12/(sqrt(2)*15)
sigma <- 15
```
```{r, echo = FALSE}
load("safeVignetteData/saviTDesignObj.RData")
```
```{r, eval=FALSE}
designObj <- designSaviT(deltaMin=deltaMin, alpha=alpha,
                         power=power, sigma=sigma,
                         alternative="greater", 
                         testType="paired", seed=1, pb=FALSE)
```
```{r}
designObj
```
The design object defines both the parameter gMom that will be used to compute the e-value, e.g. `r designObj$parameter`, and the planned sample size(s) under optional stopping, e.g. `r designObj$nPlan`. Hence, in this case we need the pre- and post-measurements of about `r designObj$nPlan[1]` patients to detect a true standarised effect size of $\delta_{\min} = 0.57$. This nPlan of `r designObj$nPlan[1]` 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 `r designObj$nPlanBatch[1]` patients (thus `r designObj$nPlanBatch[1]-designObj$nPlan[1]` more) are needed to detect $\delta_{\min} = 0.57 $ with 80% chance. <!-- in a population with $\sigma=15$. --><!--The batch planned sample size, e.g. `designObj$nPlanBatch[1]`, is non-random and serves as an upper bound for the nPlan reported above, see the [section on optional stopping](#optStopT) for further details.-->

## 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:

```{r, eval = FALSE}
designObj2 <- designSaviT(deltaMin=deltaMin, alpha=alpha,
                          nPlan=c(30, 30), sigma=sigma,
                          alternative="greater", 
                          testType="paired", seed=1, pb=FALSE)
```
```{r, echo = FALSE}
load("safeVignetteData/saviTDesignObj2.RData")
```
```{r}
designObj2
```
This reveals that due to budget constraints the experiment only has about `r round((designObj2$power)*100, 1)`% 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: 

<!-- Suppose that the tolerable type I and type II error rates are given, and that due to budget constraints we have a maximum sample size nMax of, say, 50. A prospective futility analysis can then be done by calling the design function with the provided quantities as follows: -->

```{r, eval = TRUE}
# Recall:
# alpha <- 0.05
# power <- 0.8
designObj3 <- designSaviT(nPlan=c(50, 50), 
                          alpha=alpha, power=power,
                          sigma=sigma,
                          alternative="greater",
                          testType="paired")
designObj3
```
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 `r designObj3$esMin` with 80% chance. This `r designObj3$esMin`, 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 `r designObj3$esMin`, then we can again prospective decide against conducting this experiment, or we can obtain more funds. 

# 2. <a name = "optStopT">Inference with savi tests: Optional stopping</a>
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](#optContT) 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. 

<!--In other words, this savi test retains its type I error guarantee regardless of how many times, and when the test is conducted, and irrespective of the experiment being data-dependently stopped or continued (discussed in the [next section](#optContT)). A savi test of level $\alpha = 0.05$ rejects the null whenever the observed e-value is larger than $20 = 1/\alpha$. -->

<!--TODO(Alexander): CHECK STRUCTURE 
Secondly, we show that there is a high chance of stopping early whenever the true effect size is at least as large as the minimal clinically relevant effect size.
-->

## 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. 

```{r}
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)
```
or equivalently with syntax closely resembling the standard t.test syntax:

```{r}
savi.t.test(x=preData, y=postData, 
            designObj=designObj, paired=TRUE)
```
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:

```{r, echo=FALSE}
nSim <- 1000
load("safeVignetteData/eValuesTSimple.RData")
```

```{r, eval=FALSE}
# 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}
)
```
```{r}
mean(eValues >= 20)
mean(eValues >= 20) <= alpha
```

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 `r sum(eValues > 20)` out of `r nSim` 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
<!--### Data under the alternative: <a name = "optStopT">Stopping the experiment early</a>-->
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/\alpha$, the experiment can be stopped early. This puts fewer patients at risk, and allows for more efficient scientific scrutiny. 

> Note that <a name = "martingale">not all E-variables necessarily allow for optional stopping</a>: 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](https://www.universiteitleiden.nl/binaries/content/assets/science/mi/scripties/statscience/2019-2020/thesis_rjturner_for_publication.pdf), Chapter 5.

## <a name="optionalStopPValues">Optional stopping is problematic for p-values</a>
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 `r nSim` experiments and each data set is generated with a true effect size set to zero. We first collect the t-statistics across the `r nSim` data sets, and time/sample size `r designObj$nPlan[1]`

```{r, label=zMatrix}
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$. 

```{r, label=pValuesFprCalculation}
# 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 
```
```{r, eval=FALSE}
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)
```
![](safeVignetteData/fprPValues.png)

After the second evaluation a total of `r pValueFalseRejects[3]*nSim` out of nSim=`r nSim` experiments yielded a significant p-value after the first or second observation, that is, `r pValueFalseRejects[3]*100`, 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](#optionalContinuationTypeIa). 

eValueFalseRejects
```{r, echo=FALSE}
load("safeVignetteData/eValueFalseRejectsTSimple.RData")
load("safeVignetteData/eValueFalseRejectsT.RData")
```
```{r, eval=FALSE}
# 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 
```
```{r, eval=FALSE}
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)
```
![](safeVignetteData/fprEValues.png)

Note that optional stopping always increases the chance of observing a false detection. For the savi test this increased to `r eValueFalseRejects[n1]*100`%, 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 `r pValueFalseRejects[n1]*100`%, which is well above 5%. 

```{r eval=FALSE}
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)
```
![](safeVignetteData/fprAll.png)

## Savi tests detect the effect early if it is present: deltaTrue equal to deltaMin
<!--#### Optional stopping: True effect size equals minimal clinically relevant effect size-->
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 `r nSim` 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:

```{r, echo = FALSE}
load("safeVignetteData/simDeltaTrueIsDeltaMin.RData")
```
```{r, eval = FALSE}
simDeltaTrueIsDeltaMin <- 
  sampleStoppingTimesSaviT(
    deltaTrue=deltaMin, alternative="greater", 
    testType="paired", sigma=sigma, 
    nMax=designObj$nPlan, seed=1,
    parameter=designObj$parameter,nSim=nSim)
```
```{r}
mean(
  simDeltaTrueIsDeltaMin$eValuesStopped >=
    20)
```
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:

```{r eval=FALSE}
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="")
```
![](safeVignetteData/histDeltaTrueIsDeltaMin.png)

The histogram shows the full distribution of the times at which the experiment is stopped. For instance, `r mean(simDeltaTrueIsDeltaMin$stoppingTimes < n1/2)*nSim` out of the `r nSim` 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=`r designObj$nPlan[1]`. To see the distributions of stopping times of only the experiments where the null is rejected, we run the following code:

```{r eval=FALSE}
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="")
```
![](safeVignetteData/histDeltaTrueIsDeltaMin2.png)

This can also be visualised as follows:

```{r eval=FALSE}
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)
```
![](safeVignetteData/tpr80Power.png)

Here the horizontal line represents the 80% targeted power, which is indeed reached at n=`r designObj$nPlan[1]` under optional stopping. 

## Savi tests detect the effect even earlier if deltaTrue is larger than deltaMin
<!--#### Optional stopping: True effect size larger than the minimal clinically relevant effect size-->
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

```{r, echo = FALSE}
load("safeVignetteData/simDeltaTrueLargerDeltaMin.RData")
```
```{r, eval = FALSE}
simDeltaTrueLargerDeltaMin <- 
  sampleStoppingTimesSaviT(
    deltaTrue=1.2*deltaMin, alternative="greater",
    testType="paired", sigma=sigma, 
    nMax=designObj$nPlan, seed=1,
    parameter=designObj$parameter,nSim=nSim)
```
```{r}
mean(
  simDeltaTrueLargerDeltaMin$eValuesStopped >=
    20)
```
With a larger true effect size, the power increased to `r (1-mean(simDeltaTrueLargerDeltaMin$breakVector))*100`%. 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 `r mean(simDeltaTrueIsDeltaMin$stoppingTimes)` to `r mean(simDeltaTrueLargerDeltaMin$stoppingTimes)`. This is apparent from the fact that the histogram of stopping times is now shifted to the left:

```{r eval=FALSE}
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="")
```
![](safeVignetteData/histDeltaTrueLargerThanDeltaMin.png)

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. 

```{r eval=FALSE}
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)
```
![](safeVignetteData/tprLarger.png)

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. <a name = "optContT">Optional Continuation</a>
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_](#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 \code{design} 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. 

```{r}
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: 

```{r}
pValuesBatch1 <- numeric(nSim)

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

mean(pValuesBatch1 < alpha)
```

Hence, after a first batch of data, we get `r sum(pValuesBatch1 < alpha)` incorrect null rejections out of `r nSim` experiments (`r mean(pValuesBatch1 < alpha)*100`%).

The following code continues only the non-significant `r round(mean(pValuesBatch1 > alpha)*nSim)` experiments with a second batch of data, all also generated under the null. <!--In blue the histogram of all the non-significant results (954 experiments), and in red the histogram of p-values over the combined data set. -->

```{r}
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)
```

By selectively extending the non-significant results of the first batch with a second batch of data, we got an additional `r sum(pValuesBatch2[notRejectedIndeces] < alpha)` false rejection. This brings the collective total to `r sum(pValuesBatch2 < alpha)` of false rejections out of a total of `r nSim`  studies. That is, a false positive rate of `r mean(pValuesBatch2 < alpha)*100`%, 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. 

```{r eval=FALSE}
oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
hist(pValuesBatch1, col=freqColours[1])
abline(v=0.05)
```
![](safeVignetteData/pValueHist1.png)
```{r eval=FALSE}
oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
hist(pValuesBatch2[notRejectedIndeces], col=freqColours[2])
abline(v=0.05)
```
![](safeVignetteData/pValueHist2.png)

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. 

## <a name="optionalContinuationTypeIa">i. Optional continuation by extending the experiment does not result in savi tests exceeding the tolerable $\alpha$-level</a>
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 `r n1*10` 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. 

```{r, eval=FALSE}
# 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 
```
```{r, echo = FALSE}
load("safeVignetteData/eValueFalseRejects2T.RData")
# load("safeVignetteData/allEValuesTLarge.RData")
```

```{r eval=FALSE}
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)
```
![](safeVignetteData/continueE1.png)

The original paper on e-values by [Grunwald, de Heide and Koolen](https://arxiv.org/abs/1906.07801) and [Howard, Ramdas, McAuliffe, and Sekhon]( https://doi.org/10.1214/20-AOS1991) 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. 

```{r eval=FALSE}
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)
```
![](safeVignetteData/eDrifDownwardsAsNGrows.png)

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. 

```{r, eval=FALSE}
# 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 
```

```{r, echo = FALSE}
# load("safeVignetteData/allEValuesAltT.RData")
load("safeVignetteData/eValueCorrectRejectsT.RData")
```

```{r eval=FALSE}
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)
```
![](safeVignetteData/correctReject2.png)

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=`r designObj$nPlan[1]`. 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:

```{r eval=FALSE}
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)
```
![](safeVignetteData/driftsOfAlternativeAndNull.png)

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 <a name="optionalContinuationTypeIb">meta-analysis</a>: 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](#optionalStopPValues) as original studies, and follow-up studies also with no effect. With the aforementioned original studies we carry over a false positive rate of `r mean(eStopped >= 1/alpha)*100`%. 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.

```{r, eval=FALSE}
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)
```
```{r, echo = FALSE}
#load("safeVignetteData/allEValuesAltT.RData")
load("safeVignetteData/eValueCorrectRejectsT.RData")
load("safeVignetteData/rep2.RData")
load("safeVignetteData/rep3.RData")
load("safeVignetteData/rep4.RData")
# load("safeVignetteData/repAlt.RData")
```
```{r eval=FALSE}
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)
```
![](safeVignetteData/rep2.png)
```{r}
rep2$extraRejections
```

By selectively continuing the not rejected studies we now incurred an additional `r rep2$extraRejections` false discoveries (`r mean(rep2$eOver)*100`% 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. 

```{r, eval=FALSE}
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)
```
```{r eval=FALSE}
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)
```
![](safeVignetteData/rep3.png)
```{r}
rep3$extraRejections
```
By selectively continuing the not rejected studies we now incurred no additional `r rep3$extraRejections` false discoveries retaining the cumulative `r mean(rep3$eOver)*100`% type I error rate. 

A fourth replication round with 

```{r, eval=FALSE}
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)
```
```{r eval=FALSE}
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)
```
![](safeVignetteData/rep4.png)
```{r}
rep4$extraRejections
```
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: 

```{r eval=FALSE}
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")
```
![](safeVignetteData/repEValuesDrift.png)


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. 

```{r, eval=FALSE}
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)
```
```{r, eval=FALSE}
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")
```
![](safeVignetteData/repAlt.png)
The vertical grey line is drawn at the nPlan of the original study at which we correctly reject the null with `r mean(simDeltaTrueIsDeltaMin$eValuesStopped>1/alpha)*100`% 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.  

<!-- Hence, more 

Learning from data is a continuous process, which is based on the principle that as more data accumulate the inference becomes more precise. We hope to have shown that inference based on p-values 


Optional continuation 

We believe that as more data accumulates that inference should become more precise. This 

inference is a continuous process, and that as more data accumulates one should make more precise inference. 

We hope to have highlighted the problem of optional stopping and optional 

illustrated how both optional stopping and optional continuation causes the p-values to over-reject the null. This means that both intuitive 

As such, making an observed $p< \alpha$ less informative, as it might well be due to the null hypothesis of no effect being true. This is particularly paradoxical for p-values under optional continuation, where one would expect that more data would provides more accurate inference. The inflation of the type I error implies otherwise

Robustness to optional continuation implies that intuitive actions such as 

In particular, when the null is true, then extending the experiment retains bounds the chance 

the chance of observing a low e-value remains high, 
-->
<!--
## Subconclusion
The goal the simulation studies is to illustrate the theoretical properties of savi tests developed [here](https://arxiv.org/abs/1906.07801). 

We hope to have illustrated how intuitive actions, such as monitoring the evidence and stopping early, or acquiring additional data causes the p-value to over-reject the null, thus, yielding incorrect inferences. This means that regardless of the null being true, or that there was actually an effect, the chances of 

Regardless of the presence or absence of an effect, optional stopping and optional continuation leads to an increased chance of observing $p>\alpha$, in particular, larger than $\alpha$. Especially paradoxical is the fact that optional continuation results in an over-rejection of the null, which, thus, implies that more data leads to 

As such, when confronted with $p>\alpha$ it is unclear whether the null or alternative were true. 

In contrast, when the null is true, then optional continuation actually decreases the chances of observing a large e-value, whereas if the alternative is true, then optional continuation increases the chances of observing a large e-value. This coincides with our general belief that more data provides us with more information, whereas with p-values it provides us with less information. 




the problems with classical p-values, when

the simulations show that savi tests are indeed robust to the sample sizes used in the experiment. 

that with savi tests lead to anytime valid inference, by which we mean that the conclusions remain valid regardless of one's intention to stop or to continue a study. 



Following intuition more data provides more information, whereas with p-values we get incorrect results. 

More data under the null makes chance of getting a large e-values small, whereas more data under the alternative increases the chance of observing a large e-value. 



it really doesn't matter how the 





 (n1Plan) 

on the theoretical properties 

We believe that optional stopping is the moral 

Optional continuation is of fundamental importance for science in general, as it allows us to engage to build on previous 

The latter section combined with the previous one shows e-values are robust to changes in the sample sizes of the experiment. In particular, it means that when 


This means that whenever we extend a non-rejected experiment we minimise the chance of 

This means that we reject the nullcome to the right conclusions whenever 

regardless if we 

that it does not matter whether 


monitored and 
-->

<!--
#### Subconclusion
Hence, as e-values are robust to both optional stopping and optional continuation 

Optional continuation is fundamental as it allows us to build upon previous work and 

As such, we conclude that savi tests are robust to pertubations of the experimental design. In particular, one remain convinced that 

We hope to have shed some light on how e-values are robust to the design of the experiment. This 
-->

<!-- TODO: move all two proportions instructions to separate vignette and update to the new
stream data setting
# Tests of two proportions
## 1. <a name = "design2x2">Designing Safe Experiments</a>
The safestats package also contains a safe alternative for tests of two proportions. The standard tests for this setting, which cannot deal with optional stopping, are Fisher’s exact test or the chi-squared test. These tests are applicable to data collected from two groups (indicated with "a" and "b" from here), where each data point is a binary outcome 0 (e.g. deceased) or 1 (e.g. survived). For example, group "a" might refer to the group of patients that are given the placebo, whereas group "b" is given the drug. 


### Case (I): Designing experiments with the minimal clinically relevant effect size known
As with the t-test, we might know the minimal clinically relevant effect size upfront for our test of two proportions. For example, we might only be interested in further researching or developing a drug when the difference in the proportion of cured patients in the treatment group compared to the placebo group is at least 0.3. In practice this implies, for example, that when 20% of patients get cured on average in the placebo group, we want the drug to add at least 30% to this average, so in the treated group 50% of patients should be cured. We could design a safe test for this study:

```{r }
# safeDesignProportions <- designSaviTwoProportions(deltaMin=0.3, alpha=0.05,
#                                                   power=0.80, lowN=100,
#                                                   numberForSeed = 5227)
```
For detecting this difference with a power of at least 80%, while testing at significance level 0.05, we would need:

```{r }
#safeDesignProportions$n.star
```

patients. 

A safe test could now be performed with this design object; for this, some mock data are generated below:

```{r }
# sampleExample <- as.table(matrix(c(10, safeDesignProportions[["na"]]-10, 40,
#                                    safeDesignProportions[["nb"]]-40), 
#                                  byrow=TRUE, nrow=2))
# colnames(sampleExample) <- c(0, 1)
# sampleExample
```

Performing the safe test:

```{r }
#saviTwoProportionsTest(x = sampleExample, testDesign = safeDesignProportions)
```

### Case (II): Minimal clinically relevant effect size unknown, but maximum number of samples known.
We might not have enough resources to fund our study to detect the minimal difference of 0.3. For example, we might only have funds to treat 50 patients in each group, so 100 in total. If this is the case, we could, [just as with the t-test](#powerplotT), inspect the minimal number of patients we need for the experiment to achieve a power of 80% at our significance level per effect size of interest:

```{r }
# plotResult <- plotSaviTwoProportionsSampleSizeProfile(alpha=0.05,
#                                                       power=0.80,
#                                                       highN=200, 
#                                                       maxN=100,
#                                                       numberForSeed=5222)
```

Observe that the smallest absolute difference detectable with our available resources is 0.4; we might have to cancel the study, or try to acquire more research funds, as with our current funds, we can not guarantee a high enough power for detecting the difference between the groups we are interested in. This implies that, when a non-significant result is obtained, we would be unsure whether this was caused by our underpowered study, or because there was really no difference between the groups.

Furthermore, the plot also shows the _expected sample sizes_ under optional stopping. The plot function generates experiments based on the minimal difference corresponding to the x-axis and carries out a simulation with _optional stopping_, i.e. experiments were stopped early as soon as $E >= 1/\alpha=20$ was observed, and the realised average number of patients was calculated. Observe that the difference between the planned sample size and the sample size under optional stopping is substantial. In the next section, the behaviour of the safe test for two proportions and Fisher's exact test under optional stopping is studied further.

## 2. Inference with Safe Tests and <a name = "optStop2x2">Optional Stopping</a>
#### True effect size equals minimal clinically relevant effect size
As with the safe t-test, the safe test for two proportions can be used in the optional stopping setting while retaining the type I error guarantee. In the figure below the spread of the stopping times among 1,000 simulated experiments is depicted, if the real effect size equals the minimal clinically relevant effect size as planned:

```{r }
# set.seed(5224)
# 
# optionalStoppingTrueMeanIsDesign <- 
#   simulateSpreadSampleSizeTwoProportions(
#     safeDesign=safeDesignProportions, M=1000,
#     parametersDataGeneratingDistribution=c(0.3, 0.6))
# 
# plotHistogramDistributionStoppingTimes(
#   optionalStoppingTrueMeanIsDesign, 
#   nPlan=safeDesignProportions[["n.star"]], 
#   deltaTrue = 0.3)
```

We designed the safe test such that we had a minimal power of 0.8, would the data truly come from a distribution with an absolute difference of 0.3 between the proportions of cured patients in the groups. Has this power been achieved?

```{r }
#power achieved:
#mean(optionalStoppingTrueMeanIsDesign$rejected == 1)
```

#### True effect size larger than the minimal clinically relevant effect size
We have designed the safe test for a minimal clinically relevant effect size, but what would happen if the difference between the groups was even larger in reality, i.e. if the drug had an even bigger effect?

```{r }
# set.seed(5224)
# 
# optionalStoppingTrueDifferenceBig <- 
#   simulateSpreadSampleSizeTwoProportions(
#     safeDesign=safeDesignProportions, M=1000, 
#     parametersDataGeneratingDistribution = c(0.2, 0.9))
# 
# plotHistogramDistributionStoppingTimes(
#   optionalStoppingTrueDifferenceBig, nPlan=safeDesignProportions[["n.star"]],
#   deltaTrue = 0.7)
```

We would stop, on average, even earlier! The power of the experiment also increases:

```{r }
#power achieved:
#mean(optionalStoppingTrueDifferenceBig$rejected == 1)
```

#### Data under the null: True effect size is zero, thus, much smaller than the minimal clinically relevant effect size
We can also illustrate what would happen under optional stopping, when our *null hypothesis* that there is no difference between the effect of the drug and the placebo is true:

```{r }
# set.seed(5224)
# 
# optionalStoppingTrueMeanNull <- 
#   simulateSpreadSampleSizeTwoProportions(
#     safeDesign=safeDesignProportions, M=1000, 
#     parametersDataGeneratingDistribution = c(0.5, 0.5))
# 
# plotHistogramDistributionStoppingTimes(
#   optionalStoppingTrueMeanNull, 
#   nPlan=safeDesignProportions[["n.star"]], 
#   deltaTrue = 0)
```

The type I error rate has stayed below 0.05:

```{r }
# The rate of false null rejections remained under alpha=0.05
#mean(optionalStoppingTrueMeanNull$rejected == 1)
```

#### Classical test "Fisher's exact test" under the null with optional stopping
Optional stopping, however, causes Fisher's exact test to over-reject the null. When the null is true, the rate of incorrect null rejections exceeds the tolerable $\alpha$-level:

```{r }
# set.seed(5224)
# 
# fisher_result <- simulateFisherSpreadSampleSizeOptionalStopping(
#   deltaDesign=0.5, alpha=0.05, nDesign=safeDesignProportions$n.star, 
#   power=0.8, M=100, parametersDataGeneratingDistribution=c(0.5, 0.5))
# 
# mean(fisher_result$rejected == 1)
```
Thus, 20% which is four times as much as promised. 

## 3. <a name = "optCont2x2">Optional Continuation </a>for tests of two proportions
In each of the simulations above, a fraction of the experiments did not lead to the rejection of the null hypothesis. Since safe tests allow for optional continuation, one could decide to plan a replication experiment after such a 'failed' first experiment, for example when the e-value looks promisingly high. The resulting e-values from these replication studies could then be multiplied to calculate a final e-value.

We are now going to zoom in on two of the optional stopping simulations we carried out above, where the true difference between the groups equalled our design difference (0.3), and where the true difference equalled 0. In the experiment where the true difference was 0.3, we did not reject the null in 13.2% of the studies. If we now imagine the situation we would encounter in reality, where we would not know that we were really sampling from the alternative hypothesis, how high should e-values then be to support starting a replication study? To give us some handles, we could look at the spread of e-values from studies where the null was not rejected, from our experiments under the null and under the alternative:

```{r}
# notRejectedIndex <- which(optionalStoppingTrueMeanIsDesign$rejected==FALSE)
# eValuesNotRejected <- optionalStoppingTrueMeanIsDesign$s_values[notRejectedIndex]
# nullNotRejectedIndex <- which(optionalStoppingTrueMeanNull$rejected == FALSE)
# eValuesNotRejectedNull <- optionalStoppingTrueMeanNull$s_values[nullNotRejectedIndex]
```
```{r, echo = FALSE}
# trueHist <- graphics::hist(x = eValuesNotRejected, plot = FALSE)
# nullHist <- graphics::hist(x = eValuesNotRejectedNull, plot = FALSE)
# yMax <- max(trueHist[["counts"]], nullHist[["counts"]])
# graphics::par(cex.main=1.5, mar=c(5, 6, 4, 4)+0.1, mgp=c(3.5, 1, 0), cex.lab=1.5,
#               font.lab=2, cex.axis=1.3, bty="n", las=1)
# graphics::plot(nullHist, xlim = c(0, max(eValuesNotRejected, eValuesNotRejectedNull)), 
#                freq = FALSE, col = "blue", density = 20, angle = 45, xlab = "e-values", 
#                main = "Histogram of e-values where null not rejected")
# graphics::plot(trueHist, add = TRUE, freq = FALSE, col = "red", density = 20, 
#                angle = -45)
# graphics::legend(x = "topright", legend = c("True delta: null", "True delta: design"), fill = c("blue", "red"))
```

It can be observed that, when the true difference between the groups equals our design difference, the e-values are spread out between 0 and 13. On the other hand, with our experiment under the null, all e-values were smaller than 8.

Based on this plot we could for example conclude that studies that yielded a final e-value between 10 and 20 look promising; under the null hypothesis, such high e-values were not observed in the spread plot! What would happen if we followed these studies up with a small extra study with 40 participants, and combined the resulting e-values? How many of the initially futile experiments will now lead to rejection of the null hypothesis?

```{r optionalContinuation2x2}
# continueIndex <- which(optionalStoppingTrueMeanIsDesign$s_values < 20 & 
#                          optionalStoppingTrueMeanIsDesign$s_values > 10)
# 
# interestingEValues <-
#   optionalStoppingTrueMeanIsDesign$s_values[continueIndex]
# 
# newEValues <- 
#   simulateOptionalContinuationTwoProportions(
#     interestingEValues, nFollowUp=40, 
#     parametersDataGeneratingDistribution=c(0.3, 0.6))
# 
# mean(newEValues>=20)
```

What happens when we apply this optional continuation when the data are truly generated under the null hypothesis? (note that we relax our bound of initial 'interesting' e-values here to 1, otherwise there would be no e-values to continue with)

```{r optionalContinuation2x2Null}
# continueIndex <- optionalStoppingTrueMeanNull$s_values < 20 & 
#   optionalStoppingTrueMeanNull$s_values > 1
# 
# interestingEValues <-optionalStoppingTrueMeanNull$s_values[continueIndex]
# 
# newEValues <- 
#   simulateOptionalContinuationTwoProportions(
#     interestingEValues, nFollowUp=40, 
#     parametersDataGeneratingDistribution=c(0.5, 0.5))
# 
# mean(newEValues>=20)
```
We still keep our type-I error probability guarantee.

## Short examples of usage of other testing scenarios for two proportions
Some short examples with code snippets for other testing scenarios are illustrated.

#### One-sided testing
Safe tests for two proportions can also be designed for one-sided testing. For the case when one hypothesizes that the population mean of group "a" is higher than the population mean of group "b":

```{r }
# safeDesignProportionsOneSided <- 
#   designSaviTwoProportions(deltaMin=0.5, alternative="greater",
#                            numberForSeed = 291202)
```

We can now simulate data that fit our hypothesis (more 1s observed in group "a" than in "b"):

```{r }
# sampleExampleGreater <- 
#   as.table(matrix(c(5, safeDesignProportionsOneSided[["na"]]-5, 19,
#                     safeDesignProportionsOneSided[["nb"]]-19), 
#                   byrow=TRUE, nrow=2))
# 
# colnames(sampleExampleGreater) <- c(0,1)
# sampleExampleGreater
```

This yields a high e-value:

```{r }
# saviTwoProportionsTest(x=sampleExampleGreater, 
#                        testDesign=safeDesignProportionsOneSided)
```

But if we now observe the opposite, more 1s in group "b" than in "a", the e-value will be low;

```{r }
# sampleExampleLesser <- 
#   as.table(matrix(c(safeDesignProportionsOneSided[["na"]]-5, 5,
#                     safeDesignProportionsOneSided[["nb"]]-19, 19), 
#                   byrow=TRUE, nrow=2))
# 
# colnames(sampleExampleGreater) <- colnames(sampleExampleLesser) <- c(0,1)
# sampleExampleLesser
```

```{r }
# saviTwoProportionsTest(x=sampleExampleLesser,
#                        testDesign=safeDesignProportionsOneSided)
```

#### Unbalanced design: unequal group sizes
When a balanced design is not possible, a safe test of two proportions for unequal sample sizes can be designed as well; the final ratio between the sample sizes one is going to collect has to be known for this.

```{r }
# safeDesignProportionsImbalanced <- 
#   designSaviTwoProportions(deltaMin=0.3, alpha=0.05, power=0.80, lowN=120,
#                            sampleSizeRatio=2)
# safeDesignProportionsImbalanced
```
untill here -->
