Advantages of Sequential Hypothesis Testing: 2. Flexibility and Safety

Can we early stop the testing if the p-value already reached below $\alpha$?

Again, let’s start with a simple coin toss example where we observe a sequence of independent observations $X_1,X_2,\cdots \in \{0,1\}$ from a Bernoulli distribution $B(p)$ with $p\in (0,1)$. We are interested in testing whether the underlying coin is pair ($H_0:p=0.5$) or biased toward to head ($H_1:p>0.5$).

From the previous post, we know that the binomial test is the best fixed sample size test. Formally, let’s set up the binomial test of level $\alpha =0.05$ with the minimum power of $\beta =0.95$ at $p=0.6$. In this case, the minimum sample size is given by the following R script.

set.seed(1)
printf <- function(...) invisible(print(sprintf(...)))
n_to_power <- function(n, p0, p1, alpha) {
  n <- ceiling(n)
  thres <- qbinom(alpha, n, p0, lower.tail = FALSE) # Reject null if s > thres
  pwr <- pbinom(thres, n, p1, lower.tail = FALSE)
  return(pwr)
}

power_to_n <- function(beta, p0, p1, alpha) {
  f <- function(n) {
    beta - n_to_power(n, p0, p1, alpha)
  }
  root <- stats::uniroot(f, c(1, 1e+4), tol = 1e-6)
  return(ceiling(root$root))
}

alpha <- 0.05
beta <- 0.95
p0 <- 0.5
p1 <- 0.6
n <- power_to_n(beta, p0 , p1 , alpha)
printf("%i is the minimum sample size to achieve test level %.2f and power %.2f at p= %.1f against the null p=%.2f.", 
       n, alpha, beta, p1, p0)

## [1] "280 is the minimum sample size to achieve test level 0.05 and power 0.95 at p= 0.6 against the null p=0.50."

In the standard fixed sample size setting, we compute the p-value once collecting all 280 samples. Then, if the p-value is less than the test level $\alpha$, we reject the null hypothesis. In this case, we can control the type-1 error below $\alpha$ and achieve the minimum power $\beta$ under the alternative ($p=0.6$), as the following simulation demonstrates.

binom_p_value <- function(s, n, p0) {
   return(pbinom(s-1, n, p0, lower.tail = FALSE))
}
run_binom_simul <- function(p_true, n, p0, alpha, max_iter) {
  # We shorten the simulation by using the fact X_1 + ...+X_n ~ Binom(n, p)
  s_vec <- rbinom(max_iter, size = n, prob = p_true)
  p_vec <- sapply(s_vec, binom_p_value, n = n, p0 = p0)
  return(mean(p_vec <= alpha))
}
# Under the null
max_iter <- 1e+4L
type_1_err <- run_binom_simul(p_true = p0, n, p0, alpha, max_iter)
printf("Type-1 error of the binomial test is %.2f.", type_1_err)

## [1] "Type-1 error of the binomial test is 0.04."

# Under the alternative
pwr <- run_binom_simul(p_true = p1, n, p0, alpha, max_iter)
printf("Power of the binomial test at p=%.1f is %.2f.", p1, pwr)

## [1] "Power of the binomial test at p=0.6 is 0.95."

What if we observe the coin toss one by one and we compute the p-value at each time whenever a new coin toss happens? In the previous post, we noticed that Wald’s sequential probability ratio test (SPRT) can achieve a high sample efficiency by adaptively stopping the test procedure. Why not do the same trick to the binomial test? Can we early stop the testing if the p-value already reached below $\alpha$? Let’s run a simulation to check whether this early stopping still results in a type-1 error below the test level $\alpha$.

run_binom_early_stop <- function(p_true, n, p0, alpha) {
  s <- 0
  for (i in 1:n) {
    s <- s + rbinom(1, size = 1, prob = p_true)
    p <-  binom_p_value(s, i, p0)
    if (p <= alpha) {
      return(TRUE)
    }
  }
  return(FALSE)
}
# Under the null
early_stopp_type_1_err <- mean(replicate(max_iter,{
  run_binom_early_stop(p_true = 0.5, n, p0, alpha)
  }))
printf("Type-1 error of the early stopped binomial test is %.2f.", early_stopp_type_1_err)

## [1] "Type-1 error of the early stopped binomial test is 0.26."

It turns out that the early stopping strategy increased the type 1 error significantly above the target level. This inflation of type-1 error is an example of p-hacking. This issue is getting worse if we run a longer test.

set.seed(4)
s_vec <- cumsum(rbinom(n, size = 1, prob = 0.5))
f <- function(i) binom_p_value(s_vec[i], i, p0)
p_vec <- sapply(seq_along(s_vec), f)
plot(seq_along(p_vec), p_vec, type = "l", xlab= "n", ylab ="p-value")
abline(h = alpha, col = 2)

Fig. 1 An example sample path of p-values hitting the threshold $\alpha$ although the underlying samples generated under the null. The red horizontal line corresponds to the test level $\alpha$.

Can we continue to collect more samples if the p-value at the end is slightly above $\alpha$?

Okay, we learned that we should not early stop and wait to collect all samples of the pre-calculated size. But what if the p-value is slightly above $\alpha =0.05$ - let’s say $p=0.15$. Maybe, if we would have collected a bit more samples then the p-value might possibly be less than $\alpha$ and we could reject the null. Can we continue to collect more samples to see whether the p-value becomes less than $\alpha$ or goes away from it?

Let’s run a simple simulation to answer the question. For each run, we first collect the minimum sample size $n$ computed to achieve the minimum power $\beta =0.95$. Then, compute the p-value based on the collected $n$ samples. If the p-value is less than $0.2$ then we collect another $n$ samples and re-compute the p-value based on $2n$ samples. Finally, we reject the null if the final p-value is less than or equal to the test level $\alpha$.

run_binom_a_bit_more_simul <- function(p_true, n, p0, alpha, max_iter) {
  # We shorten the simulation by using the fact X_1 + ...+X_n ~ Binom(n, p)
  s_vec <- rbinom(max_iter, size = n, prob = p_true)
  p_vec <- sapply(s_vec, binom_p_value, n = n, p0 = p0)
  
  # If p-value > alpha and < 0.2, collect n more samples
  above_alpha_ind <- which(p_vec > alpha & p_vec < 0.2)
  s_vec[above_alpha_ind] <- 
    s_vec[above_alpha_ind] + rbinom(length(above_alpha_ind), size = n, prob = p_true)
  p_vec[above_alpha_ind] <- sapply(s_vec[above_alpha_ind], binom_p_value, n = 2*n, p0 = p0)
  return(mean(p_vec <= alpha))
}
a_bit_more_type_1_err <- run_binom_a_bit_more_simul(p_true = p0, n, p0, alpha, max_iter)
printf("Type-1 error of the binomial test with the contional continuation is %.3f.", a_bit_more_type_1_err)

## [1] "Type-1 error of the binomial test with the contional continuation is 0.063."

In this case, the type-1 error increased above the test level $\alpha$ though it is not as dramatic as the early stopping case. If we continuously collect and compute p-values rather than collect a single batch then the inflation of the type-1 error becomes more significant. This type of continuation of the testing procedure is also an example of p-hacking.

run_binom_a_bit_more2_run <- function(p_true, n, p0, alpha) {
  s <- rbinom(1, size = n, prob = p_true)
  p <-  binom_p_value(s, n, p0)
  if (p <= alpha) {
    # If first n samples give p-value less than alpha
    # Reject the null
    return(TRUE)
  } else if (p >= 0.2) {
    # If first n samples give p-value above 0.2
    # Fail to reject the null and stop the test
    return(FALSE)
  }
  # Otherwise, continue the test until we collect another n samples.
  for (i in 1:n) {
    s <- s + rbinom(1, size = 1, prob = p_true)
    p <-  binom_p_value(s, i, p0)
    if (p <= alpha) {
      return(TRUE)
    }
  }
  return(FALSE)
}
# Under the null
a_bit_more2_type_1_err <- mean(replicate(max_iter,{
  run_binom_a_bit_more2_run(p_true = 0.5, n, p0, alpha)
  }))
printf("Type-1 error of the binomial test with the additional continous monitoring is %.2f.", a_bit_more2_type_1_err)

## [1] "Type-1 error of the binomial test with the additional continous monitoring is 0.18."

Always-valid p-values - flexibly and safely early stop or continue the test.

So far, we have observed how the early stop or continual test results in the type-1 error inflation for fixed sample size tests. The root cause of this inflation is that the p-value from a fixed sample size test is valid only at the pre-specified fixed sample size. In other words, let $p_n$ be the p-value of $n$ samples. Then we have $$P_{H_0}(p_n\le \alpha )\le \alpha ,\mathrm{\forall }n.$$ However, to early stop or continue the test without the type-1 error inflation, we need the following stronger inequality. $$P_{H_0}(\mathrm{\exists }n:p_n\le \alpha )\le \alpha .$$ Note that fixed sample based p-values do not necessarily satisfy the above strong inequality.¹

If a random sequence $\{p_n{\}}_{n\ge 1}$ satisfies the above stronger inequality then it is called “always-valid p-value”. Every sequential hypothesis test has the corresponding always-valid p-value. Therefore, we can flexibly and safely early stop or continue the test.

As a concrete example, let’s compute the always-valid p-value of the Wald’s SPRT. Recall that for our Bernoulli example with $H_0:p=p_0$ vs $H_1:p=p_1$, the Wald’s SPRT is a sequential hypothesis testing procedures defined as follows. * Set $S_0:=0$ and for each observation $X_n$, compute the probability ratio ${\mathrm{\Lambda }}_n$ by $${\mathrm{\Lambda }}_n:={\left(\frac{p_1}{p_0}\right)}^{X_n}{\left(\frac{1-p_1}{1-p_0}\right)}^{1-X_n}.$$

• Update $S_n:=S_{n-1}+\log {\mathrm{\Lambda }}_n$

• Make one of three following decisions:

  • If $S_n\ge b$ then stop and reject the null (since there is only one alternative, it is equivalent to accept the alternative).

  • If $S_n\le a$ then stop and reject the alternative (or accept the null).

  • Otherwise, continue to the next iteration.

Here, we can set $a=\log (1-\beta )$ and $b=\log \frac{1}{\alpha }$. In this case, we can construct an always-valid p-value by $$p_n:=min\{1,e^{-S_n}\}.$$ We can easily check $p_n\le \alpha \iff S_n\ge b=\log (1/\alpha )$. Therefore, in terms of the rejection of the null, tracking $S_n$ with the threshold $\log (1/\alpha )$ of the Wald’s SPRT is equivalent to tracking $p_n$ with the threshold $\alpha$. From the type-1 error control of the Wald’s SPRT, we can check $p_n$ is an always-valid p-value satisfying the above stronger inequality.

set.seed(1)
sprt_p_value <- function(p_true, 
                         p0 = 0.5,
                         p1 = 0.6,
                         max_iter = 1e+3L){
  s <- 0
  # Set a placeholder only for visualization.
  # We let the p-value path reaches the maximum to see the convergence
  s_history <- rep(NA, max_iter) 
  p_val_history <- rep(NA, max_iter) 
  for (n in 1:max_iter) {
    # Observe a new sample
    x <- rbinom(1, 1, p_true)
    # Update S_n
    s <- s + ifelse(x == 1, log(p1/p0), log((1-p1)/(1-p0)))
    s_history[n] <- s
    p_val_history[n] <- min(c(1,exp(-s)))
  }
  return(list(s = s_history, p_val = p_val_history))
}

# Under the null
null_out <- sprt_p_value(p_true = 0.5)
plot(seq_along(null_out$p_val), null_out$p_val, type = "l",
     xlab = "n", ylab = 'p-value',
     ylim = c(0, 1),
     main = "Under the null")
abline(h = alpha, col = 2)

# Under the alternative
alter_out <- sprt_p_value(p_true = 0.6)
plot(seq_along(alter_out$p_val), alter_out$p_val, type = "l",
     xlab = "n", ylab = 'p-value',
     ylim = c(0, 1),
     main = "Under the alternative")
abline(h = alpha, col = 2)

As a remark, the always-valid p-value above is not uniformly distributed. In fact, we can prove that $p_n\to 1$ almost surely under the null and $p_n\to 0$ almost surely under the alternative as $n\to \mathrm{\infty }$. As it is not uniformly distributed, the standard interpretation of the strength of evidence in terms of p-value is not applicable. In fact, though tracking $S_n$ and $p_n$ are equivalent, it is recommended to track $S_n$ directly as the expectation of $S_n$ corresponds to the information gain about the difference between null and alternative distributions.

a <- log(1 - beta)
b <- log(1 / alpha)
# Under the null
plot(seq_along(null_out$s), null_out$s, type = "l",
     xlab = "n", ylab = expression('S'['n']),
     main = "Under the null")
abline(h = b, col = 2)

# Under the alternative
plot(seq_along(alter_out$s), alter_out$s, type = "l",
     xlab = "n", ylab = expression('S'['n']),
     main = "Under the alternative")
abline(h = b, col = 2)