---
title: "II. Interval estimation: Classical confidence intervals and anytime-valid confidence sequences"
author: "Alexander Ly, Rosanne Turner, Udo Boehm, Rianne de Heide, Peter Grunwald"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{II. Interval estimation: Classical confidence intervals and anytime-valid confidence sequences}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

Previously we considered the testing problem. We showed that within a tolerable type I error rate alpha, we have sufficient evidence to reject the null, and stop the experiment as soon as e > 1/alpha.

Here we consider the problem of interval estimation. An anytime-valid 95\% interval estimate should contain the true value with at least 95\% chance --regardless of when, why or even if data collection stopped. 

We illustrate the following: 

1. [Classical 95\% confidence intervals fail to cover the true mean with 95\% chance across time](#classicalConfidenceIntervals)
2. [The 95\% anytime-valid confidence sequence/intervals will **always** cover the true mean with at least 95\% chance](#avConfidenceSequence)

For this demonstration we use functions from the safestats package and the following colours:

```{r, eval=FALSE}
library("safestats")
# freqColours <- c("#E31A1CE6", "#FB9A9980")
freqColours <- c("#FFDA1A", "#DAA52066")
eColours <- c("#1F78B4E6", "#A6CEE380")
alpha <- 0.05
```

# Interval estimation for a normal mean 
As a running example we consider the problem of interval estimation of a normal mean with known variance, say, $\sigma=1$. That is, we want to give an uncertainty estimate of $\mu$ where 

\begin{align}
X_{i} \overset{\text{iid}}{\sim} \mathcal{N}(\mu, \sigma^{2}).
\end{align}

<!-- Both the anytime-valid confidence interval and credible interval depend on a Gaussian/normal distribution that is centred at zero and a variance g. To highlight their difference we clarify how they act very differently with varying mixture/prior variance g. -->

For both examples we use the `computeConfidenceIntervalZ` function from the safestats package, which produces (1) classical confidence intervals if `eType="freq"`, and (2) an anytime-valid confidence sequence if `eType="mom"`. 

To assess the quality of these intervals across time, we fix a true data generating mean, say, muTrue=8, and we set sigmaTrue=1. For simplicity we take n=100 and the following code generates a single data set that we use to highlight the difference between the standard confidence intervals and the anytime-valid confidence sequence. 

```{r, eval=FALSE}
n <- 100
muTrue <- 8
sigmaTrue <- 1
  
set.seed(4)
someDataSingle <- rnorm(n, mean=muTrue,
                        sd=sigmaTrue)
```

After each observation a confidence interval can be computed yielding n=100 intervals. If any of these n=100 intervals does not contain muTrue we say that the sequence of intervals does not cover across time. 

# 1. <a name = "classicalConfidenceIntervals">Classical 95\% confidence intervals fail to cover the true mean with 95\% chance across time</a>
The classical confidence interval for a normal mean is directly related to the p-value test z-test, where the z-statistic is defined as 

\begin{align}
Z := \frac{\sqrt{n} (\bar{x} - \mu)}{\sigma}
\end{align}

All parameter values $\mu$ in a classical 95\% confidence intervals are all those $\mu$ that did not lead to a rejection at level alpha=5\%. "Inverting" the z-test shows that

\begin{align}
\bar{X} \sim \mathcal{N}(\mu, \frac{1}{n}) .
\end{align}

A 95\% confidence interval can be derived from quantiles of a normal distribution with location $\bar{x}$ and variance $1/n$.

<!-- Using the fact that $(\bar{x}-\mu)^{2}= (\mu - \bar{x})^{2}$ yields -->

### One sample path
To assess the performance of the classical confidence interval on the data generated with muTrue=8, we run the following code: 


```{r, eval=FALSE}
freqCi <- matrix(nrow=n, ncol=2)

meanVector <- 1/(1:n)*cumsum(someDataSingle)

for (i in 1:n) {
  tempResult <- computeConfidenceIntervalZ(nEff=i,
                                           meanObs=meanVector[i],
                                           eType="freq",
                                           parameter=1)
  freqCi[i, ] <- tempResult
}

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();

plot(1:n, freqCi[, 1], xlim=c(0, n), ylim=c(7, 9),
     type="l", xlab = "", ylab = "", cex.lab = 1.3, 
     cex.axis = 1.3)

polygon(c(1:n, n:1), c(freqCi[, 2], rev(freqCi[, 1])),
        col=freqColours[2], border=freqColours[1], lwd=2,
        density = NULL, angle = -20)
lines(c(1, n), c(muTrue, muTrue), lwd=2)
```
The following code shows that between the 21st and 45th observation the true mean is outside the intervals. 

```{r, eval = FALSE}
# The indices where the lower bound is above muTrue
which(freqCi[, 1] > muTrue)

# No indices where the upper bound is below muTrue
which(freqCi[, 2] < muTrue)
```

Hence, the classical confidence intervals lead to incorrect conclusions regarding muTrue, if data collection happened to stop in this range. 

### The behaviour of classical confidence intervals across time under multiple use 
The previous code worked for a specific data set/outcomes of a single experiment. Note that this 1 data set yielded 100 intervals. We now repeat the procedure, but over 1,000 data sets/experiments. The code below stores 1,000 times 100 classical confidence intervals in the variable `allFreqCis`. 
We will see that far fewer than 950 out of these 1,000 data sets/repeated experiments will cover the true mean across time. 

The following code stores a 1,000 times 100 confidence intervals: 

```{r, eval=FALSE}
mIter <- 1000

set.seed(1)
allData <- matrix(rnorm(mIter*n, mean=muTrue), nrow=mIter)
allFreqCis <- array(dim=c(mIter, n, 2))


# This indicates whether a simulation yielded an interval that 
# does not cover the true mean
freqCiError <- integer(mIter)

# This indicates the first time an interval does not cover the true mean
firstPassageTime <- rep(Inf, times=mIter)

for (sim in 1:mIter) {
  someData <- allData[sim, ]
  meanVector <- 1/(1:n)*cumsum(someData)
  
  for (i in 1:n) {
    tempResult <- computeConfidenceIntervalZ(
      nEff=i, meanObs=meanVector[i], eType="freq",
      parameter=1)
    allFreqCis[sim, i, ] <- tempResult
    
    if ((tempResult[1] > muTrue || tempResult[2] < muTrue)
        && freqCiError[sim] != 1) {
      freqCiError[sim] <- 1
      firstPassageTime[sim] <- i
    }
  }
}
```
To count the number of data sets/experiments that yielded an interval that did not cover the true mean at observation 1 to 100, we run the following code:

```{r, eval=FALSE}
complementCoverageRate <- numeric(n)

for (i in 1:n) {
  complementCoverageRate[i] <- mean(firstPassageTime <= i)
}

freqCoverageRate <- 1-complementCoverageRate

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();

plot(1:n, 100*freqCoverageRate, type="l", lwd=2, xlab="n",
     ylab="Coverage rate (%)", col=freqColours[1], ylim=c(60, 100))
lines(c(1, n), 100*c(1-alpha, 1-alpha), lwd=2, lty=2)
```
This plot reiterates the point that classical methods (i.e. p-values and confidence intervals) are made for one-shot inference. After the first time the confidence interval is computed the tolerable 5\% error was already "spent". 

The following code can be run to inspect one of the failed intervals:

```{r, eval=FALSE}
failedIndeces <- which(freqCiError==1)

someIndex <- failedIndeces[3]

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();

plot(NULL, xlim=c(0, n), ylim=c(7, 9),
     type="l", xlab = "", ylab = "", cex.lab = 1.3, 
     cex.axis = 1.3)

polygon(c(1:n, n:1), c(allFreqCis[someIndex, , 2],
                       rev(allFreqCis[someIndex, , 1])),
        col=freqColours[2], border=freqColours[1], lwd=2,
        density = NULL, angle = -20)
lines(c(1, n), c(muTrue, muTrue), lwd=2)
```

# 2. <a name = "avConfidenceSequence">The 95\% anytime-valid confidence sequence/intervals will **always** cover the true mean with at least 95\% chance</a>

The anytime-valid confidence sequence **inverts** the e-variable based on a so-called non-local moment distribution mixture. In other words, the 95\% anytime-valid confidence sequence contain all parameter values mu for which the e-variable is __less__ than 1/alpha = 20. Because the e-variable itself retains a type I error rate across time, so will the resulting confidence sequence. 

```{r, eval=FALSE}
designObj <- designSaviZ(meanDiffMin=0.5,
                         testType="oneSample",
                         sigma=sigmaTrue)
designObj
```

### One sample path
To assess the performance of the anytime-valid confidence sequence on the data with muTrue=8, we run the following code: 

```{r, eval=FALSE}
anytimeCi <- matrix(nrow=n, ncol=2)

meanVector <- 1/(1:n)*cumsum(someDataSingle)

for (i in 1:n) {
  tempResult <- computeConfidenceIntervalZ(
    nEff=i, parameter=designObj$parameter, 
    meanObs=meanVector[i])
  anytimeCi[i, ] <- tempResult
}

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();

plot(NULL, xlim=c(0, n), ylim=c(7, 9),
     type="l", xlab = "", ylab = "", cex.lab = 1.3, 
     cex.axis = 1.3)
polygon(c(1:n, n:1), c(anytimeCi[, 2], rev(anytimeCi[, 1])),
        col=eColours[2], border=eColours[1], lwd=2,
        density = NULL, angle = -20)
lines(c(1, n), c(muTrue, muTrue), lwd=2)
```
The following code shows that the anytime-valid confidence sequence covers the true mean at each moment in time. 

```{r, eval = FALSE}
# No indices where the lower bound is above muTrue
which(anytimeCi[, 1] > muTrue)

# No indices where the upper bound is below muTrue
which(anytimeCi[, 2] < muTrue)
```
The anytime-valid intervals would thus lead to the correct conclusion that muTrue is a viable value for the true mean, regardless if, whether, or why data collection is stopped. To compare the two intervals we run the following code. 

```{r, eval=FALSE}
oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();

plot(NULL, xlim=c(0, n), ylim=c(7, 9),
     type="l", xlab = "", ylab = "", cex.lab = 1.3, 
     cex.axis = 1.3)

polygon(c(1:n, n:1), c(anytimeCi[, 2], rev(anytimeCi[, 1])),
        col=eColours[2], border=eColours[1], lwd=2,
        density = NULL, angle = -20)
polygon(c(1:n, n:1), c(freqCi[, 2], rev(freqCi[, 1])),
        col=freqColours[2], border=freqColours[1], lwd=2,
        density = 60, angle = -20)
lines(c(1, n), c(muTrue, muTrue), lwd=2)
```

### The behaviour of anytime-valid confidence sequences across time under multiple use 
The previous code worked for a specific data set/outcomes of a single experiment. Note that this 1 data set yielded 100 intervals. We now repeat the procedure, but over 1,000 data sets/experiments. The code below stores 1,000 times 100 credible intervals in the variable `allAnytimeCis`. 

We will see that more than 950 out of these 1,000 data sets/repeated experiments will cover the true mean across time. 

```{r, eval=FALSE}
mIter <- 1000

set.seed(1)
allData <- matrix(rnorm(mIter*n, mean=muTrue), nrow=mIter)
allAnytimeCis <- array(dim=c(mIter, n, 2))


# This indicates whether a simulation yielded an interval that 
# does not cover the true mean
anytimeCiError <- integer(mIter)

# This indicates the first time an interval does not cover the true mean
firstPassageTime <- rep(Inf, times=mIter)

for (sim in 1:mIter) {
  someData <- allData[sim, ]
  meanVector <- 1/(1:n)*cumsum(someData)
  
  for (i in 1:n) {
    tempResult <- computeConfidenceIntervalZ(
      nEff=i, meanObs=meanVector[i],
      parameter=designObj$parameter)
    
    allAnytimeCis[sim, i, ] <- tempResult
    
    if ((tempResult[1] > muTrue || tempResult[2] < muTrue) 
        && anytimeCiError[sim] != 1) {
      anytimeCiError[sim] <- 1
      firstPassageTime[sim] <- i
    }
  }
}
```
To count the number of data sets/experiments that yielded an interval that did not cover the true mean at observation 1 to 100, we run the following code:

```{r, eval=FALSE}
complementCoverageRate <- numeric(n)

for (i in 1:n) {
  complementCoverageRate[i] <- mean(firstPassageTime <= i)
}

anytimeCiCoverageRate <- 1-complementCoverageRate

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();

plot(1:n, 100*anytimeCiCoverageRate, type="l", lwd=2,
     xlab="n", ylab="Coverage rate (%)", col=eColours[1],
     ylim=c(60, 100))
lines(c(1, n), 100*c(1-alpha, 1-alpha), lwd=2, lty=2)
```
This plot shows that the anytime-valid confidence sequences cover the true mean under repeated use as promised. 

The following code can be run to inspect one of the failed intervals:

```{r, eval=FALSE}
failedIndeces <- which(anytimeCiError==1)

someIndex <- failedIndeces[3]

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();

plot(NULL, xlim=c(0, n), ylim=c(7, 9),
     type="l", xlab = "", ylab = "", cex.lab = 1.3, 
     cex.axis = 1.3)

polygon(c(1:n, n:1), c(allAnytimeCis[someIndex, , 2],
                       rev(allAnytimeCis[someIndex, , 1])),
        col=eColours[2], border=eColours[1], lwd=2,
        density = NULL, angle = 20)
lines(c(1, n), c(muTrue, muTrue), lwd=2)
```
The following code can be used to see all coverage rate profiles. 

```{r, eval=FALSE}
oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();

plot(1:n, 100*anytimeCiCoverageRate, type="l", lwd=2,
     xlab="n", ylab="Coverage rate (%)", col=eColours[1],
     ylim=c(60, 100))
lines(1:n, 100*freqCoverageRate, col=freqColours[1], lwd=2)
lines(c(1, n), 100*c(1-alpha, 1-alpha), lwd=2, lty=2)
```

# Sensitivity wrt to the parameter gMom: Width of the confidence intervals versus anytime-valid confidence intervals

We provide some further elaboration on how/why an anytime-valid confidence sequence improves over the classical intervals. Of importance is that the widths of the classical confidence intervals is proportional to $1/\sqrt{n}$. On the other hand, the width of the anytime-valid confidence interval is stretched by the square root of a logarithmic function of $n$. To see this we plot the widths of the intervals as a function of n. 

```{r, eval=FALSE}
someA <- 0

nDomain <- 1:30
anytimeCiWidth <- freqCiWidth <- numeric(30)

for (i in nDomain) {
  tempResult <- computeConfidenceIntervalZ(
    nEff=i, meanObs=8,  eType="freq", parameter=1)
  freqCiWidth[i] <- tempResult[2]-tempResult[1]
  
  tempResult <- computeConfidenceIntervalZ(
    nEff=i, meanObs=8, parameter=designObj$parameter)
  anytimeCiWidth[i] <- tempResult[2]-tempResult[1]
}

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
plot(nDomain, anytimeCiWidth, col=eColours[1], type="l",
     lwd=2, ylim=c(0, 12))
lines(nDomain, freqCiWidth, col=freqColours[1], lwd=2,
      lty=2)
```
Hence, the anytime-valid confidence sequence is wider than the standard interval. Equivalently, the standard confidence interval is too narrow for it to be anytime-valid. 

For anytime-valid confidence sequences the following holds: For gMom small, the interval becomes infinitely wide resulting in more uncertainty. As gMom increases the width also increases, though, slowly. The following code plots the width of an anytime-valid confidence interval as a function of gMom. The following code illustrates the dependency of the interval widths on gMom. 

```{r, eval=FALSE}
gDomain <- seq(0.01, 100, by=0.01)
anytimeCiWidth <- freqIntervalWidth <- numeric(30)

for (i in seq_along(gDomain)) {
  tempResult <- computeConfidenceIntervalZ(
    nEff=1, meanObs=8, a=0, g=gDomain[i], 
    eType="freq", parameter=1)
  
  freqIntervalWidth[i] <- tempResult[2]-tempResult[1]
  
  tempResult <- computeConfidenceIntervalZ(
    nEff=1, meanObs=8, parameter=gDomain[i])
  anytimeCiWidth[i] <- tempResult[2]-tempResult[1]
}

oldPar <- setSafeStatsPlotOptionsAndReturnOldOnes();
plot(gDomain, anytimeCiWidth, col=eColours[1], type="l",
     lwd=2, ylim=c(0, 16), xlab="gMom", ylab="Interval width")
lines(gDomain, freqIntervalWidth, col=freqColours[1],
      lwd=3)
```
The plot provides an indication on how the widths of the two intervals differ for a range of hyper parameter values gMom. The increase in width of the anytime-valid confidence sequence is hardly noticeable, because it increases very slowly (logarithmically).