II. Interval estimation: Classical confidence intervals and anytime-valid confidence sequences

In safestats: Safe Anytime-Valid Inference

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
2. The 95\% anytime-valid confidence sequence/intervals will always cover the true mean with at least 95\% chance

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

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}

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.

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. Classical 95\% confidence intervals fail to cover the true mean with 95\% chance across time

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

One sample path

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

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.

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

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:

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:

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. The 95\% anytime-valid confidence sequence/intervals will always cover the true mean with at least 95\% chance

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.

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:

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.

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

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.

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:

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:

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.

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.

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.

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

safestats documentation built on Sept. 6, 2026, 1:06 a.m.