Chapter 14
Experiments, adaptive trials and decisions

Lectures 31–34 · A/B tests at scale, sequential and adaptive trials, decisions under uncertainty, and Bayesian answers with frequentist guarantees, including anytime-valid tests and confidence sequences

Chapter 13 ended with a lesson that every causal analysis eventually learns: when it can be done, randomise. A randomised experiment makes assignment ignorable by design, so the hardest assumption of observational work is simply true. This chapter is about what comes next, which is where most applied Bayesian work in industry, medicine and economics actually lives. Experiments are run in their thousands by technology companies, which call them A/B tests; they are monitored while they run, stopped early, and adapted as evidence arrives, as in modern platform trials of treatments for COVID-19; and they exist to inform decisions, which need more than a posterior: they need consequences.

The chapter also takes up an old argument. Regulators, journal editors and many statisticians judge procedures by their long-run error rates: how often a trial declares a useless drug effective. Bayesian posteriors make no such promise on their own. The last lecture shows that the two views can be reconciled, and that a Bayesian design can be tuned, by simulation or exact calculation, to have whatever frequentist guarantees are required, as the most famous vaccine trial of our time was.

• Lecture 31 is about A/B tests at scale. Over twenty-seven thousand headline experiments from the Upworthy Research Archive; a hierarchical model for the distribution of effects across experiments, with normal and heavy tails; shrinkage and the winner's curse; and how many headlines to test, following "A/B testing with fat tails".
• Lecture 32 is about sequential and adaptive trials. The peeking problem; Bayesian monitoring and predictive probabilities; adaptive allocation with Thompson sampling and its critics; and the REMAP-CAP platform trial.
• Lecture 33 is about decisions. Loss functions and Bayes actions; a trial that "failed" with $p=0.09$ and its Bayesian reanalysis; the expected value of information; and a decision in development economics.
• Lecture 34 is about frequentist guarantees for Bayesian answers. Calibrated Bayes; the Pfizer–BioNTech trial's Bayesian design and its type I error; e-values and confidence sequences, which remain valid however often you look; and Efron's frequentist accuracy of a Bayesian estimate.

Lecture 31 — A/B tests at scale

14.1 Twenty-seven thousand experiments

Between 2013 and 2015 the news site Upworthy ran an experiment for almost every story it published. Editors wrote several "packages" (a headline, sometimes with a different image) and the site showed each package to a random subset of visitors, counting impressions and clicks; the package with the most clicks became the story's headline. In 2021 J. Nathan Matias and colleagues released the complete record as the Upworthy Research Archive: 32,487 experiments, published in Scientific Data under an open licence (Matias et al., 2021). The archive was split by its creators into an exploratory sample, for developing ideas, and a larger confirmatory sample, held back so that ideas could be tested on data that did not suggest them. We follow the same discipline: we fit models on the exploratory experiments and check them on the confirmatory ones.

library(tidyverse)
library(rstan)
library(posterior)
library(jsonlite)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

harbour <- c(ink = "#10191d", muted = "#4a5f68", grid = "#d5e3e6",
             sky = "#1d5187", amber = "#805907", rose = "#ad3b28",
             moss = "#1d6b54", lilac = "#5b3a94")
theme_harbour <- function(base_size = 11) {
  theme_minimal(base_size = base_size) +
    theme(text = element_text(colour = harbour[["ink"]]),
          axis.text = element_text(colour = harbour[["muted"]]),
          panel.grid.major = element_line(colour = harbour[["grid"]]),
          panel.grid.minor = element_blank(),
          strip.text = element_text(colour = harbour[["ink"]], hjust = 0),
          legend.position = "top")
}
theme_set(theme_harbour())

Each experiment compares several packages, usually four or five. To keep the analysis simple we take two packages at random from each experiment, call them A and B, and summarise the comparison by the difference in log-odds of a click, $\widehat{d}=\text{logit}({\widehat{p}}_B)-\text{logit}({\widehat{p}}_A)$, with its usual standard error $s=\sqrt{1/c_A+1/(n_A-c_A)+1/c_B+1/(n_B-c_B)}$, where $n$ counts impressions and $c$ clicks. A difference of 0.1 in log-odds is roughly a 10% relative difference in the click rate.

osf <- c(exploratory = "3vqmp", confirmatory = "vy8mj")
dir.create("data/upworthy", showWarnings = FALSE)
for (s in names(osf)) { f <- sprintf("data/upworthy/upworthy-%s.csv", s)
  if (!file.exists(f)) download.file(paste0("https://osf.io/download/", osf[[s]], "/"), f, mode = "wb") }
pairs_of <- function(file, seed) {                       # one random pair of packages per experiment
  set.seed(seed)
  read_csv(file, show_col_types = FALSE, col_select = c(clickability_test_id, impressions, clicks)) |>
    filter(impressions >= 500, clicks > 0) |>
    group_by(clickability_test_id) |> filter(n() >= 2) |> slice_sample(n = 2) |>
    mutate(arm = c("A", "B")) |> ungroup() |>
    select(test = clickability_test_id, arm, n = impressions, c = clicks) |>
    pivot_wider(names_from = arm, values_from = c(n, c)) |>
    mutate(d = qlogis(c_B / n_B) - qlogis(c_A / n_A),
           se = sqrt(1 / c_A + 1 / (n_A - c_A) + 1 / c_B + 1 / (n_B - c_B)),
           ctr = (c_A + c_B) / (n_A + n_B))
}
ex <- pairs_of("data/upworthy/upworthy-exploratory.csv", 14)
cf <- pairs_of("data/upworthy/upworthy-confirmatory.csv", 15)
bind_rows(exploratory = ex, confirmatory = cf, .id = "sample") |>
  summarise(experiments = n(), impressions = median(c(n_A, n_B)), click_rate = median(ctr),
            median_se = median(se), significant = mean(abs(d / se) > 1.96), .by = sample) |>
  mutate(across(where(is.double), ~ round(.x, 4)))

# A tibble: 2 x 6
  sample       experiments impressions click_rate median_se significant
  <chr>              <int>       <dbl>      <dbl>     <dbl>       <dbl>
1 exploratory         4868        3125     0.013      0.219       0.242
2 confirmatory       22704        3117     0.0131     0.218       0.244

After removing packages with fewer than 500 impressions or no clicks, 4,868 exploratory and 22,704 confirmatory experiments have at least two packages to compare. A typical package was shown about 3,100 times and clicked by 1.3% of the people who saw it, so a typical comparison has a standard error of 0.22 on the log-odds scale: individual experiments are noisy. Even so, about 24% of the comparisons are "significant" at the 5% level, far more than the 5% expected if headlines made no difference. Headlines matter; the question is how much.

14.2 How large are the effects?

An observed difference ${\widehat{d}}_k$ is the true difference ${\delta }_k$ plus noise of known size, ${\widehat{d}}_k\sim \text{normal}({\delta }_k,s_k)$. With thousands of experiments, the true differences themselves have a distribution, ${\delta }_k\sim G$, and estimating $G$ is a multilevel model of the kind built in Chapter 6, with experiments as groups. Its shape matters as much as its spread. If $G$ is normal, big differences between headlines are vanishingly rare; if it has heavy tails, most headlines are about equally good but a few are far better or worse. We fit two versions: a normal $G$ with standard deviation $\sigma$, and a Student-$t$ $G$ with scale $\sigma$ and degrees of freedom $\nu$, where small $\nu$ means heavy tails.

Because only the two parameters of $G$ are of interest, we integrate out the 4,868 true differences instead of sampling them: the marginal density of each observed difference, $p({\widehat{d}}_k\mid \sigma ,\nu )=\int \text{normal}({\widehat{d}}_k\mid \delta ,s_k)G(\delta \mid \sigma ,\nu )d\delta$, is computed on a fine grid of $\delta$ values. (Sampling the latent differences directly works for the normal model, but for the heavy-tailed one HMC mixes badly, a reminder from Chapter 5 that how a model is written can matter as much as what it says.)

deconv_code <- "
data { int<lower=1> K; vector[K] d; vector<lower=0>[K] se; int<lower=1> J; vector[J] grid; int<lower=0, upper=1> heavy; }
transformed data {
  real h = grid[2] - grid[1]; matrix[K, J] kern;                    // kern[k, j] = normal(d_k | grid_j, se_k)
  for (k in 1:K) for (j in 1:J) kern[k, j] = exp(normal_lpdf(d[k] | grid[j], se[k]));
}
parameters { real<lower=0> sigma; real<lower=1, upper=100> nu; }
transformed parameters {
  vector[J] g;                                                        // G on the grid, times the grid spacing
  for (j in 1:J) g[j] = exp((heavy ? student_t_lpdf(grid[j] | nu, 0, sigma) : normal_lpdf(grid[j] | 0, sigma)) + log(h));
}
model {
  sigma ~ normal(0, 0.5); nu ~ gamma(2, 0.1);
  target += sum(log(kern * g));                                      // log of the marginal density of each d_k
}"
deconv_model <- stan_model(model_code = deconv_code)
grid <- seq(-2.5, 2.5, by = 0.02)
fit_G <- function(heavy) { path <- sprintf("fits/ch14_upworthy_%s.rds", c("normal", "heavy")[heavy + 1])
  if (file.exists(path)) return(readRDS(path))
  f <- sampling(deconv_model, data = list(K = nrow(ex), d = ex$d, se = ex$se, J = length(grid), grid = grid, heavy = heavy),
                chains = 4, iter = 1000, seed = 14, refresh = 0, pars = c("sigma", "nu"))
  saveRDS(f, path); f }
G_normal <- fit_G(0); G_heavy <- fit_G(1)
rbind(normal = summary(G_normal, pars = "sigma")$summary[, c("mean", "2.5%", "97.5%", "Rhat")],
      heavy_sigma = summary(G_heavy, pars = "sigma")$summary[, c("mean", "2.5%", "97.5%", "Rhat")],
      heavy_nu = summary(G_heavy, pars = "nu")$summary[, c("mean", "2.5%", "97.5%", "Rhat")]) |> round(3)

             mean  2.5% 97.5%  Rhat
normal      0.299 0.289 0.308 1.000
heavy_sigma 0.234 0.219 0.249 1.001
heavy_nu    4.739 3.868 5.969 1.000

Both models fit cleanly ($\widehat{R}=1.00$). The normal model puts the standard deviation of true differences at $\sigma =0.30$, so a typical pair of headlines differs by about 0.3 in log-odds, a click rate roughly a third higher for one than the other. The heavy-tailed model has a smaller scale, $\sigma =0.23$, and $\nu =4.7$ degrees of freedom (95% interval 3.9 to 6.0): most pairs of headlines are closer together than the normal model says, and a few are much further apart. A Student-$t$ with about five degrees of freedom has tails far heavier than a normal, with a variance that only just exists.

The confirmatory experiments are a fair test of the two shapes, because they played no part in fitting them. For each confirmatory comparison we compute the predictive density of the observed difference under each fitted $G$, averaging over the posterior of $\sigma$ and $\nu$, and compare the totals, exactly as LOO did in earlier chapters but with genuinely new data.

pred_logdens <- function(fit, heavy, data, draws = 100) {
  dr <- as_draws_df(fit)[round(seq(1, 2000, length.out = draws)), ]
  logk <- dnorm(outer(data$d, grid, "-"), 0, data$se, log = TRUE)
  dens <- sapply(seq_len(draws), function(i) {
    lg <- if (heavy) dt(grid / dr$sigma[i], dr$nu[i], log = TRUE) - log(dr$sigma[i]) else dnorm(grid, 0, dr$sigma[i], log = TRUE)
    m <- logk + matrix(lg + log(0.02), nrow(logk), length(grid), byrow = TRUE)
    mx <- apply(m, 1, max); mx + log(rowSums(exp(m - mx))) })
  mx <- apply(dens, 1, max); mx + log(rowMeans(exp(dens - mx)))       # log of the posterior-averaged density
}
lp_normal <- pred_logdens(G_normal, FALSE, cf); lp_heavy <- pred_logdens(G_heavy, TRUE, cf)
z <- cf$d / cf$se
c(confirmatory_experiments = nrow(cf), elpd_normal = sum(lp_normal), elpd_heavy = sum(lp_heavy),
  difference = sum(lp_heavy - lp_normal), se_difference = sd(lp_heavy - lp_normal) * sqrt(nrow(cf)),
  share_beyond_4_se_observed = mean(abs(z) > 4)) |> round(3)

  confirmatory_experiments                elpd_normal                 elpd_heavy 
                 22704.000                 -11865.508                 -11513.644 
                difference              se_difference share_beyond_4_se_observed 
                   351.864                     34.262                      0.032 

On the 22,704 confirmatory comparisons the heavy-tailed model's total log predictive density is higher by 352, with a standard error of 34: about ten standard errors, a decisive preference on data the models had never seen. The tails are visible directly, too. If true differences were normal with $\sigma =0.30$, a difference beyond four standard errors would be uncommon; 3.2% of the confirmatory comparisons have one.

sN <- mean(rstan::extract(G_normal, "sigma")$sigma); sH <- mean(rstan::extract(G_heavy, "sigma")$sigma); nH <- mean(rstan::extract(G_heavy, "nu")$nu)
tibble(delta = seq(-2, 2, by = 0.01)) |>
  mutate(normal = dnorm(delta, 0, sN), `Student-t` = dt(delta / sH, nH) / sH) |>
  pivot_longer(-delta, names_to = "shape", values_to = "density") |>
  ggplot(aes(delta, density, colour = shape)) + geom_line(linewidth = 0.9) +
  scale_y_log10(breaks = c(1e-4, 1e-3, 1e-2, 0.1, 1), labels = c("0.0001", "0.001", "0.01", "0.1", "1")) +
  scale_colour_manual(values = c(normal = harbour[["sky"]], `Student-t` = harbour[["rose"]])) +
  labs(x = "true difference between headlines (log-odds of a click)", y = "density (log scale)", colour = NULL)

14.3 Shrinkage and the winner's curse

Every experiment on the archive produced a winner, and the winner's observed advantage over the loser is a biased estimate of its true advantage. The reason is selection: among packages whose true click rates are similar, the one that wins is usually the one that was lucky, so the winning margin includes the luck. This is the winner's curse, and it has the same cure as the extreme cancer rates of Chapter 6: shrink each estimate toward what is typical, by an amount that depends on how noisy it is. With $G$ estimated from thousands of experiments, the posterior mean of each ${\delta }_k$ is

and its shape depends on $G$. A normal $G$ shrinks every estimate by the same proportion. A heavy-tailed $G$ shrinks small estimates hard, since they are probably noise around a typical effect near zero, and leaves large estimates almost alone, since heavy tails make large true effects believable. The widget shows both.

post_mean <- function(d, se, heavy) {                          # E[delta | d] on the grid, at the posterior means of G
  lg <- if (heavy) dt(grid / sH, nH, log = TRUE) - log(sH) else dnorm(grid, 0, sN, log = TRUE)
  w <- exp(dnorm(outer(d, grid, "-"), 0, se, log = TRUE) + matrix(lg, length(d), length(grid), byrow = TRUE))
  as.vector(w %*% grid) / rowSums(w) }
winners <- ex |> mutate(winner_adv = abs(d), sign = sign(d), z = abs(d) / se) |> filter(z > 1.96)
winners <- winners |> mutate(shrunk_normal = sign * post_mean(d, se, FALSE), shrunk_heavy = sign * post_mean(d, se, TRUE))
winners |> summarise(winners = n(), observed_advantage = mean(winner_adv), shrunk_normal = mean(shrunk_normal),
                     shrunk_heavy = mean(shrunk_heavy)) |> mutate(across(-winners, ~ round(.x, 3)))
writeLines(paste0("var C14U=", toJSON(list(sN = round(sN, 4), sH = round(sH, 4), nH = round(nH, 3), ctr = round(median(ex$ctr), 4)),
                                      auto_unbox = TRUE), ";"), "data/ch14_upworthy.js")

# A tibble: 1 x 4
  winners observed_advantage shrunk_normal shrunk_heavy
    <int>              <dbl>         <dbl>        <dbl>
1    1178              0.665         0.393        0.404

Of the 4,868 exploratory comparisons, 1,178 produced a "significant" winner. Their average observed advantage is 0.67 in log-odds, about a 95% higher click rate. Under either estimate of $G$ the posterior mean advantage of the same winners is about 0.40, so roughly 40% of the observed advantage is the winner's curse. The two shapes agree on the average because most winners are of moderate size; they disagree about individual experiments, as the widget shows: under the Student-$t$ $G$ a winner that won by a large margin keeps most of it, while one that only just won is shrunk heavily.

14.4 How many headlines should we test?

The shape of $G$ is not an academic detail; it changes how experiments should be run. Suppose an editor has a fixed number of visitors to spend on testing a story's headline, writes $m$ candidate headlines, shows each to an equal share of the visitors, and then uses the one with the highest posterior mean for everyone else. More headlines give more chances of finding an exceptional one; fewer headlines give each a larger sample, so the winner is picked more reliably. The best $m$ balances the two, and the balance depends on the tails of $G$.

set.seed(14); p <- median(ex$ctr)
best_of <- function(m, visitors, draw, R = 40000) {
  se <- sqrt(1 / (visitors / m * p * (1 - p)))                 # standard error of one headline's log-odds
  D <- matrix(draw(R * m), R, m); O <- D + matrix(rnorm(R * m, 0, se), R, m)
  mean(D[cbind(1:R, max.col(O))]) }                            # true value of the headline that looked best
draws <- list(normal = function(k) rnorm(k, 0, sN), `Student-t` = function(k) sH * rt(k, nH))
lean <- expand_grid(visitors = c(5000, 20000), m = c(2, 3, 4, 6, 8, 12, 16, 24, 32), G = names(draws)) |>
  mutate(gain = pmap_dbl(list(m, visitors, G), function(m, v, g) best_of(m, v, draws[[g]])))
lean |> group_by(visitors, G) |> slice_max(gain, n = 1) |> mutate(gain = round(gain, 3)) |> print()
ggplot(lean, aes(m, gain, colour = G)) + geom_line(linewidth = 0.9) + geom_point(size = 1.4) +
  facet_wrap(~ visitors, labeller = labeller(visitors = function(v) paste(format(as.numeric(v), big.mark = ","), "visitors"))) +
  scale_x_log10(breaks = c(2, 4, 8, 16, 32)) +
  scale_colour_manual(values = c(normal = harbour[["sky"]], `Student-t` = harbour[["rose"]])) +
  labs(x = "number of headlines tested (log scale)", y = "expected gain (log-odds)", colour = NULL)

# A tibble: 4 x 4
# Groups:   visitors, G [4]
  visitors     m G          gain
     <dbl> <dbl> <chr>     <dbl>
1     5000    16 Student-t 0.304
2     5000    12 normal    0.274
3    20000    32 Student-t 0.489
4    20000    16 normal    0.405

With 5,000 visitors, the normal $G$ says that 12 headlines is best, the heavy-tailed $G$ that 16 is; with 20,000 visitors the normal optimum is 16 and the heavy-tailed one 32 or more, the largest number simulated. Under heavy tails the expected gain keeps rising with the number of headlines, because the occasional exceptional headline is so much better than the rest that it is worth trying many candidates, each on a modest sample, to find it. Under a normal $G$ the gains from more candidates level off sooner and the noise of small samples wins. Upworthy's editors tested four or five headlines per story; on these estimates, they could have gained by writing more.

Lecture 32 — Sequential and adaptive trials

14.5 The peeking problem

An A/B test on a website accumulates data continuously, and the dashboard shows the current result at every moment. It is very tempting to stop as soon as the result looks decisive. A classical significance test is not designed for this. It controls the false-positive rate at 5% for one look at a pre-specified sample size; look twenty times and stop at the first $|z|>1.96$, and a useless change will be declared a winner far more often than 5% of the time. The same is true of a naive Bayesian rule that stops when the posterior probability that B is better exceeds 0.975. The simulation below runs 20,000 tests in which the two versions are identical, checking after each of $L$ equal batches of visitors.

set.seed(14)
run_mean <- function(inc) {                                    # running mean of each row's batch estimates
  cs <- inc; if (ncol(inc) > 1) for (j in 2:ncol(inc)) cs[, j] <- cs[, j - 1] + inc[, j]
  cs / matrix(seq_len(ncol(inc)), nrow(inc), ncol(inc), byrow = TRUE) }
peek_once <- function(looks, delta, R = 20000, tau = sN, batch_se = 0.15) {  # one batch alone has standard error batch_se
  inc <- matrix(rnorm(R * looks, rep(delta, each = 1), batch_se), R, looks)       # per-batch estimates of the difference
  est <- run_mean(inc)                                                            # running mean
  v <- batch_se^2 / matrix(1:looks, R, looks, byrow = TRUE)                       # its variance
  z_hit <- rowSums(abs(est / sqrt(v)) > 1.96) > 0                                 # classical test at every look
  post_sd <- sqrt(1 / (1 / v + 1 / tau^2)); post_m <- post_sd^2 * est / v         # normal prior N(0, tau^2)
  pp <- pnorm(post_m / post_sd)
  b_hit <- rowSums(pp > 0.975 | pp < 0.025) > 0                                    # Bayesian rule at every look
  c(z_test = mean(z_hit), bayes_rule = mean(b_hit)) }
peeking <- map_dfr(c(1, 2, 5, 10, 20, 50), ~ c(looks = .x, peek_once(.x, delta = 0)))
peeking |> mutate(across(-looks, ~ round(.x, 3)))

# A tibble: 6 x 3
  looks z_test bayes_rule
  <dbl>  <dbl>      <dbl>
1     1  0.051      0.028
2     2  0.082      0.054
3     5  0.14       0.108
4    10  0.198      0.168
5    20  0.248      0.218
6    50  0.318      0.292

With a single look, the classical test has its promised 5% false-positive rate. With 5 looks it is 14%, with 20 looks 25% and with 50 looks 32%. The Bayesian rule, which stops when the posterior probability (under a normal prior with the standard deviation estimated from the Upworthy experiments) leaves the interval from 0.025 to 0.975, does a little better because the prior shrinks early, noisy estimates, but it too reaches 22% with 20 looks and 29% with 50. Stopping at the first exciting moment makes exciting moments easy to find.

Does this mean that Bayesian inference is ruined by optional stopping? No, but the guarantee it offers is different. The posterior is computed from the likelihood, and the likelihood does not depend on why the experimenter stopped, so a posterior probability of 0.975 means what it says averaged over the prior: among all experiments whose effects really come from the prior and that stopped with a posterior probability of 0.975 that B is better, B is better in about 97.5% of them, however the stopping was decided. What it does not promise is a 2.5% error rate for every fixed effect, including an effect of exactly zero. The next simulation draws each experiment's true difference from the normal prior and checks the claims made by the Bayesian rule with twenty looks.

set.seed(15); R <- 40000; looks <- 20; batch_se <- 0.15
delta <- rnorm(R, 0, sN)
inc <- matrix(rnorm(R * looks, rep(delta, looks), batch_se), R, looks)
est <- run_mean(inc); v <- batch_se^2 / matrix(1:looks, R, looks, byrow = TRUE)
post_sd <- sqrt(1 / (1 / v + 1 / sN^2)); pp <- pnorm(post_sd^2 * est / v / post_sd)
stop_at <- apply(pp > 0.975 | pp < 0.025, 1, function(h) if (any(h)) which(h)[1] else NA)
called <- !is.na(stop_at); pp_stop <- pp[cbind(which(called), stop_at[called])]
said_B <- pp_stop > 0.5
c(share_stopped = mean(called), mean_look_at_stop = mean(stop_at, na.rm = TRUE),
  correct_when_said_B = mean(delta[called][said_B] > 0), correct_when_said_A = mean(delta[called][!said_B] < 0)) |> round(3)

      share_stopped   mean_look_at_stop correct_when_said_B correct_when_said_A 
              0.871               3.631               0.990               0.990 

With twenty looks, 87% of the experiments stopped, on average at the fourth look. When the rule stopped and said B was better, B really was better 99.0% of the time; when it said A was better, A really was better 99.0% of the time. Both figures exceed the nominal 97.5%, because the rule stops as soon as the threshold is crossed and the posterior often overshoots it. Optional stopping did not damage the calibration of the posterior at all.

The two guarantees answer different questions. The Bayesian guarantee is the one a decision-maker wants if the prior really describes the population of ideas being tested, which is why the hierarchical model of Lecture 31 matters: the Upworthy archive tells us what that population looks like. The frequentist guarantee is the one a regulator wants, because it holds for every possible effect, including the one a sceptic believes, that the drug does nothing. Lecture 34 shows how to have both.

14.6 Monitoring a trial: posterior and predictive probabilities

Clinical trials are monitored too, but formally: a data monitoring committee looks at pre-specified interim analyses and stops the trial early if the evidence is overwhelming, if the treatment is futile, or if it is harmful. Two Bayesian quantities guide such decisions. The posterior probability of benefit, $Pr({\theta }_T<{\theta }_C\mid \text{data so far})$ for a mortality outcome, says what we believe now. The predictive probability of success says how likely the trial is to succeed if it runs to completion: average, over the posterior, the probability that the data still to come will push the final analysis over its success threshold. It is the natural quantity for futility decisions, because a trial whose predictive probability of success is 2% is not worth continuing, whatever the current estimate.

Consider an interim analysis of a two-arm trial in critically ill patients, planned for 400 patients per arm with a final success criterion of a posterior probability above 0.99 that mortality is lower on treatment. At the interim, 200 patients per arm have been followed: 60 deaths on treatment and 72 on control.

set.seed(14)
interim <- c(n_t = 200, d_t = 60, n_c = 200, d_c = 72); final_n <- 400; success <- 0.99
post_t <- function(R) rbeta(R, 1 + interim[["d_t"]], 1 + interim[["n_t"]] - interim[["d_t"]])
post_c <- function(R) rbeta(R, 1 + interim[["d_c"]], 1 + interim[["n_c"]] - interim[["d_c"]])
th_t <- post_t(4000); th_c <- post_c(4000)
p_benefit_now <- mean(th_t < th_c)
final_prob <- function(dt, dc) {                                  # posterior Pr(benefit) at the final analysis
  a <- rbeta(2000, 1 + dt, 1 + final_n - dt); b <- rbeta(2000, 1 + dc, 1 + final_n - dc); mean(a < b) }
ppos <- mean(map_lgl(1:1000, function(i) {                       # simulate the remaining 200 patients per arm
  dt <- interim[["d_t"]] + rbinom(1, final_n - interim[["n_t"]], th_t[i])
  dc <- interim[["d_c"]] + rbinom(1, final_n - interim[["n_c"]], th_c[i])
  final_prob(dt, dc) > success }))
c(mortality_treatment = 60 / 200, mortality_control = 72 / 200, posterior_prob_benefit = p_benefit_now,
  predictive_prob_success = ppos) |> round(3)

    mortality_treatment       mortality_control  posterior_prob_benefit 
                  0.300                   0.360                   0.895 
predictive_prob_success 
                  0.292 

At the interim, mortality is 30% on treatment and 36% on control, and the posterior probability that treatment lowers mortality is 0.90, which sounds encouraging. But the predictive probability of meeting the trial's demanding success criterion at the end is only 0.29: a benefit of the observed size, if real, would probably not be estimated precisely enough with 400 patients per arm to reach a posterior probability of 0.99. The two numbers answer different questions, and a monitoring committee needs both.

14.7 Adaptive allocation

If the evidence so far favours one arm, why keep giving half the patients, or half the visitors, the arm that looks worse? Response-adaptive randomisation shifts the allocation toward arms that are doing better. Its most elegant form is Thompson sampling (Thompson, 1933): for each new unit, draw one value of every arm's success rate from its current posterior, and assign the unit to the arm whose draw is largest. An arm is therefore chosen with exactly its posterior probability of being best. Early on every arm gets a fair share; as evidence accumulates, the allocation concentrates on the winner. Thompson sampling is a multi-armed bandit strategy, designed to maximise the total reward while learning, and it is how many companies now run headline and recommendation tests.

We compare it with equal allocation on simulated headline tests with four candidates, whose true differences in log-odds are drawn from the heavy-tailed distribution of Section 14.2, a base click rate of 1.4% and 20,000 visitors.

run_test <- function(true_lo, visitors = 20000, thompson = TRUE, batch = 100, drift = 0) {
  K <- length(true_lo); clicks <- shown <- numeric(K); share <- matrix(NA, visitors / batch, K)
  for (b in seq_len(visitors / batch)) {
    base <- drift * (b - 1) / (visitors / batch)                         # optional drift in the base rate over time
    if (thompson) { draws <- matrix(rbeta(K * 200, 1 + clicks, 1 + shown - clicks), 200, K, byrow = TRUE)
      w <- tabulate(max.col(draws), K) / 200 } else w <- rep(1 / K, K)
    n_b <- as.vector(rmultinom(1, batch, w)); share[b, ] <- n_b / batch
    clicks <- clicks + rbinom(K, n_b, plogis(qlogis(0.014) + true_lo + base)); shown <- shown + n_b }
  list(clicks = sum(clicks), best_found = which.max((clicks + 1) / (shown + 2)) == which.max(true_lo),
       est_gap = qlogis((clicks[which.max(true_lo)] + 1) / (shown[which.max(true_lo)] + 2)) - qlogis((sum(clicks[-which.max(true_lo)]) + 1) / (sum(shown[-which.max(true_lo)]) + 2)),
       true_gap = max(true_lo) - log(mean(exp(true_lo[-which.max(true_lo)]))), share = share)
}
set.seed(14)
sims <- map_dfr(1:300, function(i) { lo <- sH * rt(4, nH)
  e <- run_test(lo, thompson = FALSE); t <- run_test(lo, thompson = TRUE)
  tibble(design = c("equal", "Thompson"), clicks = c(e$clicks, t$clicks), best_found = c(e$best_found, t$best_found),
         gap_error = c(e$est_gap - e$true_gap, t$est_gap - t$true_gap)) })
sims |> summarise(mean_clicks = mean(clicks), best_found = mean(best_found), mean_gap_error = mean(gap_error), .by = design) |>
  mutate(across(where(is.double), ~ round(.x, 3)))
set.seed(3); one <- run_test(c(0.25, 0, -0.1, 0.05))
as_tibble(one$share, .name_repair = ~ paste("headline", 1:4)) |> mutate(visitors = row_number() * 100) |>
  pivot_longer(-visitors, names_to = "headline", values_to = "share") |>
  ggplot(aes(visitors, share, fill = headline)) + geom_area(position = "stack", alpha = 0.85) +
  scale_fill_manual(values = unname(harbour[c("rose", "sky", "moss", "amber")])) +
  labs(x = "visitors", y = "share of each batch", fill = NULL)

# A tibble: 2 x 4
  design   mean_clicks best_found mean_gap_error
  <chr>          <dbl>      <dbl>          <dbl>
1 equal           290.      0.757          0    
2 Thompson        363.      0.793         -0.078

Over 300 simulated tests, Thompson sampling earned 363 clicks per test against 290 for equal allocation, 25% more, because it stopped wasting traffic on poor headlines early. It also found the truly best headline slightly more often (79% against 76%). The price shows in the estimates: Thompson sampling underestimated the gap between the best headline and the others by 0.08 on average, while equal allocation was unbiased. Arms that fall behind early are starved of data, so their early bad luck is never corrected, and the comparison inherits it.

Adaptive randomisation is controversial in medicine, for reasons the simulation hints at (Robertson et al., 2023). Patients enrolled late are allocated differently from patients enrolled early, so if the patient population or standards of care drift during the trial, a time trend becomes confounded with treatment, and the comparison is biased unless the analysis models time. Imbalanced arms estimate the difference less precisely than equal ones, so a trial may need more patients for the same power. And the benefit to patients inside the trial is often small compared with the benefit to the far larger number treated after it, who need a precise answer. Platform trials such as REMAP-CAP use adaptive allocation with care: allocation is updated at scheduled analyses rather than patient by patient, a fixed share is kept for control, and time is modelled.

14.8 Platform trials

A platform trial studies several treatments for one disease under a single master protocol that runs indefinitely: arms are added as new candidates appear, dropped when a statistical trigger declares them effective, futile or harmful, and analysed together with a shared control group. The design was developed for exactly the conditions of a pandemic, and Bayesian models are its natural language, because they make posterior probabilities available at any moment and borrow strength across related arms.

remap <- tibble(arm = factor(c("control", "tocilizumab", "sarilumab"), c("control", "tocilizumab", "sarilumab")),
                n = c(397, 350, 45), deaths = c(142, 98, 10))
remap_code <- "
data { int<lower=1> A; array[A] int<lower=0> n; array[A] int<lower=0> deaths; }
parameters { real alpha; vector[A - 1] beta; }                 // log-odds of death, and log odds ratios against control
model { alpha ~ normal(0, 1.5); beta ~ normal(0, 1);
  deaths[1] ~ binomial_logit(n[1], alpha); deaths[2:A] ~ binomial_logit(n[2:A], alpha + beta); }"
remap_fit <- stan(model_code = remap_code, data = list(A = 3, n = remap$n, deaths = remap$deaths),
                  chains = 4, iter = 2000, seed = 14, refresh = 0)
b <- rstan::extract(remap_fit, "beta")$beta
or_survival <- exp(-b)                                           # odds ratio for survival, as REMAP-CAP reports
rbind(tocilizumab = c(median = median(or_survival[, 1]), quantile(or_survival[, 1], c(0.025, 0.975)), p_benefit = mean(b[, 1] < 0)),
      sarilumab = c(median = median(or_survival[, 2]), quantile(or_survival[, 2], c(0.025, 0.975)), p_benefit = mean(b[, 2] < 0))) |> round(3)
c(p_sarilumab_better_than_tocilizumab = mean(b[, 2] < b[, 1])) |> round(3)

            median  2.5% 97.5% p_benefit
tocilizumab  1.422 1.056 1.919     0.993
sarilumab    1.860 0.974 3.740     0.969
p_sarilumab_better_than_tocilizumab 
                              0.782 

From the published counts, with a simple model and no adjustment, the posterior median odds ratio for survival is 1.42 for tocilizumab (95% interval 1.06 to 1.92), with a posterior probability of benefit of 0.993, and 1.86 for sarilumab (0.97 to 3.74), with probability 0.969. The trial's own adjusted estimates for survival were 1.64 (1.14 to 2.35) and 2.01 (1.18 to 4.71). Ours are smaller and less certain, as expected from a model that ignores prognostic covariates (age, illness severity, site and time) which the trial adjusted for, and which sharpen comparisons in a randomised trial. The posterior probability that sarilumab is better than tocilizumab is 0.78: with only 48 patients on sarilumab, the data cannot separate the two drugs, which is why REMAP-CAP's hierarchical model pooled them as a class.

Lecture 33 — Decisions

14.9 From posteriors to actions

A posterior distribution summarises what we know; a decision also needs to know what we want. Bayesian decision theory joins the two. There is a set of possible actions $a$ (approve the drug or not, ship headline B or keep A, fund the programme or not), an unknown state $\theta$, and a loss $L(a,\theta )$, or equivalently a utility $U(a,\theta )=-L(a,\theta )$, saying how bad each action is in each state.

Familiar summaries are Bayes actions for particular losses. If the action is a number and the loss is squared error, $L(a,\theta )=(a-\theta {)}^2$, the Bayes action is the posterior mean; with absolute error it is the posterior median. The most useful case in practice is a choice between two actions with asymmetric consequences. Suppose treating a patient who would benefit gains $b$, and treating one who would not costs $c$ (side effects, money, a forgone alternative). Then treating has posterior expected utility $bPr(\text{benefit}\mid \text{data})-c(1-Pr(\text{benefit}\mid \text{data}))$, and it is the better action exactly when

The threshold is set by the costs, not by convention. A cheap, safe treatment for a deadly disease should be used at a posterior probability well below 0.5; an expensive one with serious side effects needs much more. The 0.95 and 0.975 of significance testing are not decision thresholds at all; they are conventions for controlling a particular error rate, which may or may not correspond to anybody's losses.

14.10 A trial that "failed"

The EOLIA trial randomised 249 patients with very severe acute respiratory distress syndrome to early extracorporeal membrane oxygenation (ECMO), in which a machine oxygenates the blood outside the body, or to conventional ventilation, with ECMO allowed as a rescue (Combes et al., 2018). At 60 days, 44 of 124 patients on early ECMO had died (35%), against 57 of 125 on conventional treatment (46%). The relative risk was 0.76 with a 95% confidence interval from 0.55 to 1.04 and $p=0.09$, and the trial, which had been stopped early for futility, was reported as negative.

We reproduce the reanalysis with three priors on the log relative risk, each combined with a uniform prior on the control group's risk: a minimally informative prior (normal with standard deviation 10); a sceptical prior centred on no effect that gives only a 5% chance to a relative risk below 0.67; and an enthusiastic prior centred on 0.78 that gives a 15% chance to harm.

eolia <- list(d1 = 44, n1 = 124, d0 = 57, n0 = 125)
rr_code <- "
data { int d1; int n1; int d0; int n0; real m; real s; }
parameters { real<lower=0, upper=1> p0; real<upper=-log(p0)> log_rr; }     // p1 = p0 * RR must stay below 1
transformed parameters { real p1 = p0 * exp(log_rr); }
model { log_rr ~ normal(m, s); d0 ~ binomial(n0, p0); d1 ~ binomial(n1, p1); }
generated quantities { real rr = exp(log_rr); real arr = p0 - p1; }"
rr_model <- stan_model(model_code = rr_code)
priors <- tribble(~prior, ~m, ~s,
  "minimally informative", 0, 10,
  "sceptical", 0, log(1 / 0.67) / qnorm(0.95),
  "enthusiastic", log(0.78), -log(0.78) / qnorm(0.85))
eolia_post <- priors |> mutate(draws = map2(m, s, function(m, s) {
  f <- sampling(rr_model, data = c(eolia, m = m, s = s), chains = 4, iter = 2000, seed = 14, refresh = 0)
  as_tibble(rstan::extract(f, c("rr", "arr"))) }))
eolia_post |> mutate(summ = map(draws, ~ tibble(rr_median = median(.x$rr), lo = quantile(.x$rr, 0.025), hi = quantile(.x$rr, 0.975),
    p_benefit = mean(.x$rr < 1), p_large = mean(.x$rr < 0.67), p_arr_2pct = mean(.x$arr > 0.02)))) |>
  select(prior, prior_sd = s, summ) |> unnest(summ) |> mutate(across(where(is.double), ~ round(.x, 3)))

# A tibble: 3 x 8
  prior                prior_sd rr_median    lo    hi p_benefit p_large p_arr_2pct
  <chr>                   <dbl>     <dbl> <dbl> <dbl>     <dbl>   <dbl>      <dbl>
1 minimally informati~   10         0.776 0.559  1.03     0.958   0.184      0.918
2 sceptical               0.243     0.833 0.641  1.08     0.914   0.05       0.841
3 enthusiastic            0.24      0.776 0.592  0.99     0.979   0.138      0.946

Our reanalysis reproduces the paper's main result. Under the minimally informative prior the posterior median relative risk is 0.78 (95% interval 0.56 to 1.03), the probability of any benefit is 0.96, and the probability of an absolute reduction in mortality of at least 2 percentage points is 0.92. A sceptical prior pulls the median to 0.83 and the probability of benefit to 0.91; an enthusiastic one gives 0.78 and 0.98. The probability of a large benefit, a relative risk below 0.67, is only 0.05 to 0.18 under the three priors. A clinician who reads "$p=0.09$, not significant" as "ECMO does not work" has misread a trial whose data make a moderate benefit highly probable.

eolia_post |> select(prior, draws) |> unnest(draws) |>
  ggplot(aes(rr, colour = prior)) + geom_density(linewidth = 0.9, adjust = 1.2) +
  geom_vline(xintercept = 1, linetype = 2, colour = harbour[["muted"]]) +
  scale_colour_manual(values = c(`minimally informative` = harbour[["ink"]], sceptical = harbour[["sky"]], enthusiastic = harbour[["rose"]])) +
  coord_cartesian(xlim = c(0.4, 1.4)) + labs(x = "relative risk of death, ECMO against conventional care", y = NULL, colour = NULL) +
  theme(axis.text.y = element_blank())

14.11 The value of information

Should a hospital adopt early ECMO after EOLIA, and should anyone run another trial? Decision theory answers both with the same posterior. To make the example concrete we need values that the trial does not supply, and we state them as illustrative assumptions, which a real analysis would take from health economics: each death averted gains 10 quality-adjusted life years (QALYs), a health system is willing to pay $\lambda$ per QALY, and early ECMO costs an extra $40,000 per patient. The net monetary benefit per patient of adopting ECMO is then

where $p_0-p_1$ is the absolute reduction in the risk of death. The Bayes decision adopts ECMO if the posterior mean of NB is positive. But the decision is made under uncertainty, and resolving that uncertainty has a value. The expected value of perfect information (EVPI) is how much better, on average, the decision would be if we knew $\theta$:

Multiplied by the number of patients the decision affects, EVPI is an upper bound on what any further research is worth. The expected value of sample information (EVSI) does the same calculation for a specific new trial of $n$ patients per arm: simulate the trial's data from the current posterior, update, decide, and average. If EVSI exceeds the trial's cost, the trial is worth running.

post <- eolia_post$draws[[1]]                                       # minimally informative posterior
nb <- function(arr, wtp) wtp * 10 * arr - 40000
evpi <- function(arr, wtp) mean(pmax(nb(arr, wtp), 0)) - max(mean(nb(arr, wtp)), 0)
evsi <- function(n, wtp, R = 4000) {                                 # preposterior analysis with a normal approximation
  set.seed(14); i <- sample(nrow(post), R); p0 <- 57 / 125
  arr_true <- post$arr[i]; m0 <- mean(post$arr); v0 <- var(post$arr)
  se_new <- sqrt(2 * p0 * (1 - p0) / n)                              # standard error of a new trial's risk difference
  obs <- rnorm(R, arr_true, se_new); v1 <- 1 / (1 / v0 + 1 / se_new^2); m1 <- v1 * (m0 / v0 + obs / se_new^2)
  mean(pmax(nb(m1, wtp), 0)) - max(nb(m0, wtp), 0) }
tibble(wtp = c(20000, 50000, 100000)) |>
  mutate(mean_NB = map_dbl(wtp, ~ mean(nb(post$arr, .x))), P_NB_positive = map_dbl(wtp, ~ mean(nb(post$arr, .x) > 0)),
         EVPI = map_dbl(wtp, ~ evpi(post$arr, .x)), EVSI_200 = map_dbl(wtp, ~ evsi(200, .x)),
         EVSI_1000 = map_dbl(wtp, ~ evsi(1000, .x))) |> mutate(across(-c(wtp, P_NB_positive), round), P_NB_positive = round(P_NB_positive, 3)) |>
  mutate(wtp = format(wtp, big.mark = ",", scientific = FALSE)) |> as.data.frame()
writeLines(paste0("var C14V=", toJSON(list(m = round(mean(post$arr), 5), s = round(sd(post$arr), 5)), auto_unbox = TRUE), ";"), "data/ch14_voi.js")

      wtp mean_NB P_NB_positive EVPI EVSI_200 EVSI_1000
1  20,000  -19223         0.059  339       82       243
2  50,000   11941         0.645 7058     4598      6397
3 100,000   63883         0.856 4357     1736      3591

The decision depends on how much a QALY is worth. At $20,000 per QALY the expected net benefit is negative (−$19,200 per patient) and ECMO should not be adopted; at $50,000 it is positive ($11,900) but uncertain, with a 64% chance that adoption is the right decision; at $100,000 it is clearly positive. The expected value of perfect information peaks where the decision is closest, at $7,060 per patient at $50,000 per QALY. A new trial with 200 patients per arm would recover about $4,580 of it per patient, and one with 1,000 per arm about $6,340. With thousands of patients with severe ARDS treated each year, even a fraction of that per patient amounts to tens of millions of dollars, far more than a trial costs: decision theory says that a confirmatory trial is worth running, a conclusion that "not significant" alone could not reach.

14.12 A decision in development economics

Economists face the same structure when they recommend policies. Section 6.5 analyses Rachael Meager's (2019) Bayesian hierarchical model of seven randomised evaluations of microcredit, in which loans were made available in some places and not others. Suppose a government is considering an expansion of microcredit in a new country. The relevant evidence is not the average effect across the studied sites but the posterior predictive effect in a new site, which includes the variation between sites. We take the five studies that measured household consumption, estimate each site's effect, fit the hierarchical model and compare the two distributions.

mc_url <- "https://raw.githubusercontent.com/wwiecek/baggr/master/data/microcredit_simplified.rda"
if (!file.exists("data/microcredit/microcredit_simplified.rda")) {
  dir.create("data/microcredit", showWarnings = FALSE); download.file(mc_url, "data/microcredit/microcredit_simplified.rda", mode = "wb") }
load("data/microcredit/microcredit_simplified.rda")                 # consumption, USD PPP per two weeks
sites <- c("Mexico", "Mongolia", "Bosnia", "India", "Morocco")
site_est <- microcredit_simplified |> group_by(group) |>
  summarise(tau = coef(lm(consumption ~ treatment))[2], se = sqrt(vcov(lm(consumption ~ treatment))[2, 2]), n = n()) |>
  mutate(site = sites[group])
hier_code <- "
data { int<lower=1> K; vector[K] y; vector<lower=0>[K] se; }
parameters { real mu; real<lower=0> tau; vector[K] z; }
transformed parameters { vector[K] theta = mu + tau * z; }
model { mu ~ normal(0, 50); tau ~ normal(0, 25); z ~ std_normal(); y ~ normal(theta, se); }
generated quantities { real theta_new = normal_rng(mu, tau); }"
hier_model <- stan_model(model_code = hier_code)
micro_fit <- sampling(hier_model, data = list(K = nrow(site_est), y = site_est$tau, se = site_est$se),
                      chains = 4, iter = 4000, seed = 14, refresh = 0, control = list(adapt_delta = 0.999))
md <- as_tibble(rstan::extract(micro_fit, c("mu", "tau", "theta_new")))
site_est |> select(site, n, tau, se) |> mutate(across(c(tau, se), ~ round(.x, 2)))
cost_per_year <- 50                                                  # illustrative cost per household, USD PPP per year
rbind(mu = c(mean = mean(md$mu), quantile(md$mu, c(0.05, 0.95)), p_positive = mean(md$mu > 0)),
      new_site = c(mean = mean(md$theta_new), quantile(md$theta_new, c(0.05, 0.95)), p_positive = mean(md$theta_new > 0)),
      net_benefit_new_site = c(mean = mean(26 * md$theta_new - cost_per_year), quantile(26 * md$theta_new - cost_per_year, c(0.05, 0.95)),
                               p_positive = mean(26 * md$theta_new > cost_per_year))) |> round(2)
c(divergences = get_num_divergent(micro_fit), tau_median = median(md$tau)) |> round(2)

# A tibble: 5 x 4
  site         n   tau    se
  <chr>    <int> <dbl> <dbl>
1 Mexico   16560  5.51  4.1 
2 Mongolia   961 50.4  32.2 
3 Bosnia     995 -5.17 19.2 
4 India     6827  4.64  5.48
5 Morocco   5487 -2.94  6.02
                      mean      5%    95% p_positive
mu                    3.73   -4.81  12.66       0.79
new_site              3.88  -12.04  20.43       0.72
net_benefit_new_site 50.79 -362.99 481.28       0.61
divergences  tau_median 
          0           5 

The five site estimates differ a great deal in precision: Mexico's, from 16,560 households, has a standard error of 4.1 dollars per fortnight, Mongolia's, from 961, one of 32.2. The posterior mean of the average effect is $3.73 per fortnight (90% interval −$4.81 to $12.66), with a 79% probability of being positive, close to Meager's published 3.4 for the same outcome. For a new site the predictive distribution is much wider, from −$12.04 to $20.43, because the sites differ (the posterior median of the between-site standard deviation $\tau$ is $5). With an illustrative programme cost of $50 per household per year and 26 fortnights a year, the expected net benefit in a new country is $51 per household per year, but the probability that the benefit exceeds the cost is only 0.61. A decision-maker who looked only at the average effect, which is positive with probability 0.79, would overstate the case for expansion; the honest summary is that it is a reasonable bet with a high chance of losing money.

Lecture 34 — Bayesian answers with frequentist guarantees

14.13 Calibrated Bayes

Statisticians have argued for a century about whether probability describes beliefs or long-run frequencies. In applied work the argument has largely been settled by a division of labour that Roderick Little (2006) called calibrated Bayes, building on George Box (1980) and Donald Rubin (1984): use Bayesian models to make inferences, because they handle complex models, prior information, nuisance parameters and missing data coherently; and use frequentist calculations to check them, because a procedure that is badly calibrated in repeated use should not be trusted, whatever its philosophy.

Two different guarantees are in play, and Lecture 32 has already separated them.

• Calibration on average over the prior. If parameters really are drawn from the prior and data from the model, then 90% credible intervals contain the truth 90% of the time and posterior probabilities are accurate, by construction. This is what simulation-based calibration checks in Chapter 5, where it is used to test the computation.
• Calibration for every parameter value. A frequentist guarantee, such as a type I error of at most 2.5%, must hold for each fixed $\theta$, including values the prior thinks unlikely. Bayesian procedures do not have it automatically, but it can be checked, and a design can be tuned to have it, by choosing the decision thresholds, the prior or the sample size.

The rest of this lecture shows four ways of getting the second guarantee: tuning the thresholds of a Bayesian design (Section 14.14), using Bayes factors that are valid under optional stopping (Section 14.15), turning those Bayes factors into intervals that are valid at every moment (Section 14.16), and measuring the frequentist variability of a Bayesian estimate directly (Section 14.17). Two further ideas deserve a mention. Probability-matching priors are priors chosen so that one-sided credible intervals have exact or approximate frequentist coverage; for a normal mean the flat prior is one, and for many models Jeffreys' prior comes close. Conformal Bayesian computation (Fong and Holmes, 2021) uses the posterior predictive distribution to build prediction intervals that have finite-sample frequentist coverage whatever the model, a guarantee that is especially useful when a model is surely misspecified.

14.14 Tuning a Bayesian design: the Pfizer–BioNTech trial

Chapter 1 analysed the final result of the BNT162b2 vaccine trial. Its design deserves a closer look, because it is the most widely seen example of a Bayesian trial tuned to frequentist requirements (Polack et al., 2020). The analysis concerned $\theta$, the probability that a COVID-19 case occurred in the vaccine group rather than the placebo group. With equal-sized groups, $\theta =(1-\text{VE})/(2-\text{VE})$, where VE is the vaccine efficacy, so that $\text{VE}=(1-2\theta )/(1-\theta )$. The prior was beta(0.700102, 1), chosen so that its mean corresponds to an efficacy of 30% and it is only weakly informative. The trial would declare success if the posterior probability that VE exceeded 30% was above 0.995 at any of four interim analyses, after 32, 62, 92 and 120 cases, or above 0.986 at the final analysis after 164 cases.

Why these numbers? The regulators required that the probability of declaring success when the true efficacy is only 30%, the type I error, be at most 2.5% across all five looks. A posterior probability threshold of 0.95, applied at every look, would not achieve that. The thresholds 0.995 and 0.986 were chosen to make it so. Because the number of vaccine cases at each look is binomial, the operating characteristics can be computed exactly.

a0 <- 0.700102; b0 <- 1; theta0 <- (1 - 0.3) / (2 - 0.3)           # theta at VE = 30%
looks <- c(32, 62, 92, 120, 164)
boundary <- function(thr) map2_dbl(looks, thr, function(N, t) {    # most vaccine cases still counting as success
  x <- 0:N; ok <- pbeta(theta0, a0 + x, b0 + N - x) > t; if (any(ok)) max(x[ok]) else -1 })
p_success <- function(theta, bnd) {                                 # exact probability of success at some look
  p <- 1; prev <- 0; succ <- 0
  for (j in seq_along(looks)) {
    q <- convolve(p, rev(dbinom(0:(looks[j] - prev), looks[j] - prev, theta)), type = "open")
    q <- pmax(q, 0)
    if (bnd[j] >= 0) { succ <- succ + sum(q[1:(bnd[j] + 1)]); q[1:(bnd[j] + 1)] <- 0 }
    p <- q; prev <- looks[j] }
  succ }
designs <- list(protocol = c(0.995, 0.995, 0.995, 0.995, 0.986), `0.95 at every look` = rep(0.95, 5),
                `0.95 at the final look only` = c(1, 1, 1, 1, 0.95))
bounds <- map(designs, boundary)
bounds$protocol |> setNames(paste(looks, "cases")) |> print()
ve_to_theta <- function(ve) (1 - ve) / (2 - ve)
map_dfr(designs, function(d) { b <- boundary(d)
  c(type_I_error = p_success(theta0, b), power_VE_50 = p_success(ve_to_theta(0.5), b), power_VE_60 = p_success(ve_to_theta(0.6), b)) }, .id = "design") |>
  mutate(across(-design, ~ round(.x, 4)))
th <- rbeta(1e6, a0 + 8, b0 + 162); ve <- (1 - 2 * th) / (1 - th)   # the final data: 8 vaccine cases of 170
c(post_prob_VE_above_30 = pbeta(theta0, a0 + 8, b0 + 162), VE_median = median(ve), quantile(ve, c(0.025, 0.975))) |> round(4)
writeLines(paste0("var C14P=", toJSON(list(a0 = a0, looks = looks), auto_unbox = TRUE), ";"), "data/ch14_pfizer.js")

 32 cases  62 cases  92 cases 120 cases 164 cases 
        6        15        25        35        53 
# A tibble: 3 x 4
  design                      type_I_error power_VE_50 power_VE_60
  <chr>                              <dbl>       <dbl>       <dbl>
1 protocol                          0.0218       0.458       0.884
2 0.95 at every look                0.136        0.762       0.976
3 0.95 at the final look only       0.0548       0.683       0.965
post_prob_VE_above_30             VE_median                  2.5% 
               1.0000                0.9486                0.9035 
                97.5% 
               0.9763 

The protocol's boundaries allow success with at most 6 vaccine cases among the first 32, 15 among 62, 25 among 92, 35 among 120 and 53 among 164. With these thresholds the probability of declaring success when the true efficacy is exactly 30% is 2.18%, within the required 2.5%, while the power is 46% if the efficacy is 50% and 88% if it is 60%. Using 0.95 at every look would have given a type I error of 13.6%, more than five times the limit; using 0.95 only at the final analysis would have given 5.5%, still too much because a posterior probability of 0.95 is, with this nearly flat prior, roughly a one-sided test at 5%, not 2.5%. The final data, 8 cases in the vaccine group and 162 in the placebo group, give a posterior probability above 0.9999 that efficacy exceeds 30%, and a posterior median efficacy of 94.9% with a 95% credible interval from 90.3% to 97.6%, matching the published 95.0% (90.3% to 97.6%).

This is calibration by design. The prior and the posterior probabilities make the trial's statements interpretable ("given the data, the probability that efficacy exceeds 30% is above 99.9%"), while the thresholds give regulators the error rate they require. The same recipe, specify the Bayesian decision rule, compute or simulate its operating characteristics over a range of true effects, and tune until they are acceptable, is the backbone of the US Food and Drug Administration's draft guidance on Bayesian methods in drug trials (FDA, 2026). When exact calculation is impossible, as in REMAP-CAP, the operating characteristics are found by simulating thousands of trials.

14.15 Looking as often as you like: e-values

The thresholds of Section 14.14 work for a fixed schedule of five looks. A website wants to look continuously. For that, the right object is a statistic that keeps its meaning however often it is inspected, and a Bayes factor provides one. Let ${\widehat{d}}_t$ be the estimated difference after $t$ batches, with variance $v_t$. Compare the null hypothesis $\delta =0$ with an alternative that spreads $\delta$ as normal$(0,{\tau }^2)$:

Why does this work? Recall from Chapter 1 that a martingale is a sequence whose expected next value, given everything so far, equals its current value. A Bayes factor monitored over time is one, under the null hypothesis, and a classical inequality bounds how often such a sequence can ever climb high.

For the e-value above, $E_t$ is exactly such a Bayes factor, so rejecting the null hypothesis the first time $E_t$ reaches $1/\alpha$ has a type I error of at most $\alpha$, however many times, and at whatever moments, the data are inspected, and even if the experimenter decides to continue or stop by looking at the data.

Statistics with this property are called e-values (or, as a process over time, e-processes), and inference built on them is always valid (Johari et al., 2022; Grünwald, de Heide and Koolen, 2024). The mixture Bayes factor above is the "mixture sequential probability ratio test" used by several commercial experimentation platforms. The price is power: to be valid at every moment, the boundary must be higher than a single-look test's. The simulation compares a single test at the end with an e-value monitored after every batch, under no effect and under effects drawn from the Upworthy distribution.

set.seed(14); R <- 20000; looks <- 50; batch_se <- 0.3; tau_e <- sN
run_e <- function(delta) {
  inc <- matrix(rnorm(R * looks, rep(delta, looks), batch_se), R, looks)
  est <- run_mean(inc); v <- batch_se^2 / matrix(1:looks, R, looks, byrow = TRUE)
  E <- exp(dnorm(est, 0, sqrt(v + tau_e^2), log = TRUE) - dnorm(est, 0, sqrt(v), log = TRUE))
  hit <- E >= 20; first <- apply(hit, 1, function(h) if (any(h)) which(h)[1] else NA)
  c(fixed_test_rejects = mean(abs(est[, looks] / sqrt(v[, looks])) > 1.96), e_value_rejects = mean(!is.na(first)),
    e_value_mean_looks_when_rejecting = mean(first, na.rm = TRUE)) }
rbind(no_effect = run_e(rep(0, R)), effects_from_G = run_e(sH * rt(R, nH))) |> round(3)

               fixed_test_rejects e_value_rejects
no_effect                   0.049           0.019
effects_from_G              0.738           0.632
               e_value_mean_looks_when_rejecting
no_effect                                  15.39
effects_from_G                             14.34

Under no effect, the single test at the end rejects 4.9% of the time, as designed, and the e-value monitored at all 50 looks rejects only 1.9% of the time: Ville's bound of 5% holds with room to spare, however the looks are used. When effects are drawn from the Upworthy distribution, the e-value detects them 63% of the time against 74% for the single test, but when it does reject it does so after 14 looks on average instead of 50: experiments with real effects can be stopped far earlier. That trade, a little power for the freedom to look and stop at any time, is why continuous-monitoring platforms use e-values.

The e-value above used a normal approximation to the difference in click rates. For counts there is an exact version, and it is worth seeing because it shows exactly where the guarantee comes from. Suppose two headlines are shown with fixed probabilities, a share $q_0$ of visitors seeing B. If the headlines are equally good, each click, whenever it arrives, comes from B with probability $q_0$, whatever the click rate is. So we can watch only the clicks: after $n$ clicks, $k$ of them on B. Under the null hypothesis the clicks are Bernoulli($q_0$); under the alternative the share of clicks on B is some other $q$, about which we put a beta($a$, $a$) prior. The Bayes factor after $n$ clicks is

with $B$ the beta function. The relative lift of B over A is $\frac{q}{1-q}·\frac{1-q_0}{q_0}$, so a prior on $q$ is a prior on the lift.

Two things about Ville's inequality deserve emphasis. The guarantee holds for every prior: the beta($a$, $a$) prior can be chosen to reflect the effects seen in past experiments, and that changes how quickly real effects are detected, never the error rate. And the guarantee survives any stopping rule, including stopping because the result looks good, which is exactly what experimenters do. In the language of this lecture: the prior tunes power; the martingale guarantees validity.

The Upworthy deconvolution of Section 14.2 gives a prior. Its heavy-tailed effect distribution has scale 0.234 on the log-odds scale, and a beta($a$, $a$) prior on the click share has a log-odds standard deviation of about $\sqrt{2/a}$, so $a=37$ matches it. We compare it with a uniform prior ($a=1$), first by simulation: 20,000 visitors per arm, a click rate of 1.5%, and either no effect or a lift of 20%. The Bayes factor is checked after every click and the test stops when it reaches 20 ($\alpha =0.05$); for comparison, a $z$-test on the click share is applied either once at the end or after every 5% of the clicks.

log_bf <- function(x, q0, a) { k <- cumsum(x); n <- seq_along(x)     # log Bayes factor after every click (x = 1 for B)
  lbeta(a + k, a + n - k) - lbeta(a, a) - (k * log(q0) + (n - k) * log(1 - q0)) }
one_test <- function(p_A, p_B, visitors, a) {
  x <- sample(c(rep(0, rbinom(1, visitors, p_A)), rep(1, rbinom(1, visitors, p_B))))   # clicks in random order
  n <- length(x); stop_at <- which(log_bf(x, 0.5, a) >= log(20))[1]
  looks <- unique(pmax(1, round(n * (1:20) / 20))); z <- (cumsum(x)[looks] - looks / 2) / sqrt(looks / 4)
  c(bf_rejects = !is.na(stop_at), share_of_clicks_used = ifelse(is.na(stop_at), 1, stop_at / n),
    z_once_rejects = abs(tail(z, 1)) > 1.96, z_20_looks_rejects = any(abs(z) > 1.96)) }
set.seed(3)
sims <- expand_grid(prior_a = c(1, 37), lift = c(0, 0.2)) |>
  mutate(res = map2(prior_a, lift, ~ replicate(ifelse(.y == 0, 4000, 2000), one_test(0.015, 0.015 * (1 + .y), 20000, .x)))) |>
  mutate(bf_rejects = map_dbl(res, ~ mean(.x["bf_rejects", ])),
         median_share_used_when_rejecting = map_dbl(res, ~ median(.x["share_of_clicks_used", .x["bf_rejects", ] == 1])),
         z_once_rejects = map_dbl(res, ~ mean(.x["z_once_rejects", ])),
         z_20_looks_rejects = map_dbl(res, ~ mean(.x["z_20_looks_rejects", ]))) |> select(-res)
sims |> mutate(across(where(is.double), ~ round(.x, 3))) |> as.data.frame()

  prior_a lift bf_rejects median_share_used_when_rejecting z_once_rejects
1       1  0.0      0.034                            0.043          0.047
2       1  0.2      0.242                            0.495          0.632
3      37  0.0      0.018                            0.470          0.049
4      37  0.2      0.380                            0.588          0.649
  z_20_looks_rejects
1              0.253
2              0.779
3              0.238
4              0.768

With no effect, the $z$-test checked 20 times rejects about a quarter of the time, as Section 14.5 found. The Bayes factor, checked after every one of roughly 600 clicks, rejects 3.4% of the time with the uniform prior and 1.8% with the Upworthy prior: both within Ville's 5%, as the theorem promises. With a real 20% lift, the informed prior detects the effect in 38% of experiments against 24% for the uniform prior, and when it does so it has used, at the median, 59% of the clicks. A single $z$-test at the end has more power (about 64%), and that is the trade: the fixed test must wait for all 20,000 visitors and may be looked at only once; the Bayes factor can be watched continuously and stops early when the evidence is strong, at the cost of power at the planned sample size.

Now apply the test to the 22,704 confirmatory Upworthy experiments of Section 14.1, taking one pair of headlines from each. The archive records only totals, so we put each experiment's clicks in a random order, which is how they would arrive if the null hypothesis were true and is a reasonable approximation otherwise, and use the share of impressions shown version B as $q_0$.

set.seed(14)
up <- read_csv("data/upworthy/upworthy-confirmatory.csv", show_col_types = FALSE,
               col_select = c(clickability_test_id, impressions, clicks)) |>
  filter(impressions >= 500, clicks > 0) |> group_by(clickability_test_id) |> filter(n() >= 2) |> slice_sample(n = 2) |>
  mutate(arm = c("A", "B")) |> ungroup() |> select(test = clickability_test_id, arm, n = impressions, c = clicks) |>
  pivot_wider(names_from = arm, values_from = c(n, c)) |> mutate(q0 = n_B / (n_A + n_B))
seq_tests <- pmap_dfr(up, function(test, n_A, n_B, c_A, c_B, q0) {
  x <- sample(c(rep(0, c_A), rep(1, c_B))); lb <- log_bf(x, q0, 37); stop_at <- which(lb >= log(20))[1]
  n <- c_A + c_B; z <- (c_B - n * q0) / sqrt(n * q0 * (1 - q0))
  tibble(clicks = n, bf_stops = !is.na(stop_at), share_used = ifelse(is.na(stop_at), 1, stop_at / n), z_final_significant = abs(z) > 1.96) })
tibble(experiments = nrow(seq_tests), median_clicks = median(seq_tests$clicks),
  bf_stops = mean(seq_tests$bf_stops), z_final_significant = mean(seq_tests$z_final_significant),
  bf_stops_but_z_not = mean(seq_tests$bf_stops & !seq_tests$z_final_significant),
  median_share_used_when_stopping = median(seq_tests$share_used[seq_tests$bf_stops])) |>
  pivot_longer(everything()) |> mutate(value = format(round(value, 3), big.mark = ",", drop0trailing = TRUE, scientific = FALSE)) |> as.data.frame()

                             name     value
1                     experiments    22,704
2                   median_clicks        88
3                        bf_stops     0.076
4             z_final_significant     0.249
5              bf_stops_but_z_not     0.001
6 median_share_used_when_stopping     0.664

A $z$-test at the end of each experiment finds a significant difference in 24.9% of them, which with the heavy-tailed effects of Section 14.2 is plausible: many headline pairs really do differ. The sequential Bayes factor stops 7.6% of the experiments, almost all of them among those the $z$-test also flags, and uses a median of 66% of the clicks when it stops. Most Upworthy experiments are small (the median has 88 clicks in all), and with so little data an always-valid test can rarely be sure; the fixed test, looked at once, finds three times as many differences. For a newsroom testing headlines for a few hours, the fixed test with a fixed stopping time is the better tool. For a platform that runs experiments for weeks and wants to stop as soon as a winner is clear, or to stop a harmful change early, the always-valid test is the safe one.

14.16 Confidence sequences

A test can be turned into an interval by collecting the values it does not reject. Applying the Bayes factor to every possible click share $q$, not just $q_0$, gives the set

and Ville's inequality, applied at the true $q$, says that $Pr(q\in C_n\text{ for every }n)\ge 1-\alpha$. This is a confidence sequence: a sequence of intervals that contain the truth simultaneously at all times with probability at least $1-\alpha$, so it can be reported after every click without any correction. Note the shape of the condition. The ratio is the prior predictive probability of the data divided by their likelihood at $q$, and that equals the prior density at $q$ divided by the posterior density at $q$. So $C_n$ is the set where the posterior density exceeds $\alpha$ times the prior density: a posterior region, with its level set by the prior rather than by the posterior mass. Robbins (1970) found this construction; Howard, Ramdas, McAuliffe and Sekhon (2021) and Waudby-Smith and Ramdas (2024) made it practical and sharp.

The figure follows one large Upworthy experiment click by click: the 95% confidence sequence for the relative lift, and a naive 95% interval recomputed after every click, which at the last click is the ordinary fixed-sample interval.

big <- up |> filter(c_A + c_B > 600, abs(q0 - 0.5) < 0.03) |>                # a large, evenly split experiment
  mutate(z = (c_B - (c_A + c_B) * q0) / sqrt((c_A + c_B) * q0 * (1 - q0))) |> arrange(desc(abs(z))) |> slice(1)
big |> select(n_A, n_B, c_A, c_B, z) |> mutate(z = round(z, 2)) |> as.data.frame()
set.seed(21)
x <- sample(c(rep(0, big$c_A), rep(1, big$c_B))); if (big$c_B < big$c_A) x <- 1 - x   # B is the winner
q0 <- if (big$c_B < big$c_A) 1 - big$q0 else big$q0
n <- seq_along(x); k <- cumsum(x); qg <- seq(0.02, 0.98, by = 0.0005)
lift <- function(q) (q / (1 - q)) * ((1 - q0) / q0) - 1
cs <- map_dfr(unique(round(exp(seq(log(20), log(length(x)), length.out = 120)))), function(t) {
  ok <- lbeta(37 + k[t], 37 + t - k[t]) - lbeta(37, 37) - (k[t] * log(qg) + (t - k[t]) * log(1 - qg)) < log(20)
  qhat <- k[t] / t; se <- sqrt(qhat * (1 - qhat) / t)
  tibble(t = t, cs_lo = lift(min(qg[ok])), cs_hi = lift(max(qg[ok])),
         naive_lo = lift(qhat - 1.96 * se), naive_hi = lift(qhat + 1.96 * se)) })
ggplot(cs, aes(t)) + geom_ribbon(aes(ymin = cs_lo, ymax = cs_hi), fill = harbour[["sky"]], alpha = 0.25) +
  geom_line(aes(y = naive_lo), linetype = 2, colour = harbour[["rose"]]) + geom_line(aes(y = naive_hi), linetype = 2, colour = harbour[["rose"]]) +
  geom_hline(yintercept = 0, linetype = 3) + scale_x_log10() + scale_y_continuous(labels = scales::percent) +
  coord_cartesian(ylim = c(-0.6, 1.2)) + labs(x = "clicks so far (log scale)", y = "relative lift of B over A")
tail(cs, 1) |> mutate(across(-t, ~ round(.x, 3))) |> as.data.frame()
first_excl <- c(confidence_sequence = min(cs$t[cs$cs_lo > 0]), naive = min(cs$t[cs$naive_lo > 0])); first_excl

   n_A  n_B c_A c_B     z
1 6294 6176 459 239 -8.08
    t cs_lo cs_hi naive_lo naive_hi
1 698 0.389 1.587    0.617    1.213
confidence_sequence               naive 
                162                  89 

The experiment compared two headlines shown to about 6,200 visitors each; the better one received 459 clicks against 239. After all 698 clicks, the 95% confidence sequence puts the relative lift between 39% and 159%, against 62% to 121% for the fixed-sample interval: the price of validity at every moment is an interval about twice as wide. The naive interval first excludes zero after about 90 clicks, the confidence sequence after about 160. The naive interval's early exclusion is legitimate only if 90 clicks was the sample size fixed in advance; an experimenter who watches the naive interval and stops when it first excludes zero is running the peeking procedure of Section 14.5, whose error rate is several times the nominal one. The confidence sequence can be watched without that worry.

Confidence sequences now exist for far more than proportions: for means of bounded outcomes without any distributional assumption, by "betting" (Waudby-Smith and Ramdas, 2024), for quantiles, for regression coefficients, and, using asymptotic arguments, for average treatment effects estimated with propensity scores and outcome models (Waudby-Smith, Arbour, Sinha, Kennedy and Ramdas, 2024). When the null hypothesis is composite, so that no single likelihood ratio is available, universal inference (Wasserman, Ramdas and Balakrishnan, 2020) builds e-values from a likelihood ratio computed on a split of the data, and the theory of e-processes extends Ville's inequality to these cases.

14.17 The frequentist accuracy of a Bayesian estimate

A posterior standard deviation describes uncertainty given the data. How much would the posterior mean itself vary if the study were repeated? Bradley Efron (2015) gave a general answer that needs nothing beyond the posterior draws already in hand. Write the data as $x$ and the likelihood as $f(x\mid \theta )$, and let ${\alpha }_x(\theta )={\nabla }_x\log f(x\mid \theta )$ be the gradient of the log likelihood with respect to the data. Then the gradient of the posterior expectation of any quantity $t(\theta )$ with respect to the data is a posterior covariance,

which can be estimated from posterior draws, and the delta method turns it into a frequentist standard deviation, ${\text{sd}}_{\text{freq}}\approx \sqrt{{\widehat{g}}^{\top }V\widehat{g}}$, with $\widehat{g}$ the estimated gradient and $V$ the sampling covariance of $x$. For the microcredit model of Section 14.12, the data are the five site estimates, $x_k\sim \text{normal}({\theta }_k,s_k)$, so ${\alpha }_x(\theta {)}_k=({\theta }_k-x_k)/s_k^2$ and $V=\text{diag}(s_k^2)$. We compare Efron's estimate with the posterior standard deviation, and check it with a parametric bootstrap that refits the model 200 times to data simulated from the fitted site effects.

th_draws <- rstan::extract(micro_fit, "theta")$theta; x <- site_est$tau; s <- site_est$se
alpha <- sweep(sweep(th_draws, 2, x), 2, s^2, "/")                   # (theta_k - x_k) / s_k^2 for every draw
g <- as.vector(cov(md$mu, alpha))                                    # Efron's gradient of E[mu | x]
efron_sd <- sqrt(sum(g^2 * s^2))
set.seed(14); theta_hat <- colMeans(th_draws)
boot <- map_dbl(1:200, function(i) {
  xb <- rnorm(length(x), theta_hat, s)
  fb <- sampling(hier_model, data = list(K = length(x), y = xb, se = s), chains = 2, iter = 2000, seed = i,
                 refresh = 0, control = list(adapt_delta = 0.99))
  mean(rstan::extract(fb, "mu")$mu) })
c(posterior_mean_mu = mean(md$mu), posterior_sd_mu = sd(md$mu), efron_frequentist_sd = efron_sd,
  bootstrap_sd = sd(boot)) |> round(2)

   posterior_mean_mu      posterior_sd_mu efron_frequentist_sd 
                3.73                 5.94                 2.99 
        bootstrap_sd 
                3.40 

The posterior standard deviation of $\mu$ is $5.94, but the frequentist standard deviation of its posterior mean is only about $3: Efron's formula gives $2.99 and the bootstrap $3.40. They answer different questions. The posterior standard deviation describes what we do not know about $\mu$, including what the prior and the small number of sites leave open about the variation between sites. The frequentist standard deviation describes how much the estimate would move if the five studies were repeated with the same true site effects: not much, because shrinkage stabilises it. A Bayesian estimate can be more stable than its posterior uncertainty suggests, and Efron's formula measures the stability without refitting the model; the bootstrap, which needed 200 refits, agrees reasonably well.

14.18 Experiments and decisions in practice

The ideas of this chapter come together in a few practical rules. Design experiments with their decisions in mind: state the action that each result would lead to, and the costs of each kind of error. Use what is known about the population of effects, which for organisations that run many experiments is a hierarchical model of their own history. Decide in advance how often you will look, and use a method that is valid for that schedule. Report posterior probabilities of effects of meaningful size, not only of any effect. And, where a frequentist guarantee is required, compute the operating characteristics of the Bayesian procedure and tune it until they are met.

14.19 Exercises

Each exercise has a hint and an answer button. Code for the answers is in R/ex14.R.

Further reading and references

1. Azevedo, E. M., Deng, A., Montiel Olea, J. L., Rao, J. and Weyl, E. G. (2020). A/B testing with fat tails. Journal of Political Economy, 128(12), 4614–4672.
2. Berry, S. M., Carlin, B. P., Lee, J. J. and Müller, P. (2010). Bayesian Adaptive Methods for Clinical Trials. CRC Press.
3. Box, G. E. P. (1980). Sampling and Bayes' inference in scientific modelling and robustness. Journal of the Royal Statistical Society: Series A, 143(4), 383–430.
4. Combes, A., Hajage, D., Capellier, G., et al. (2018). Extracorporeal membrane oxygenation for severe acute respiratory distress syndrome. New England Journal of Medicine, 378(21), 1965–1975.
5. Dimairo, M., Pallmann, P., Wason, J., et al. (2020). The Adaptive designs CONSORT Extension (ACE) statement: a checklist with explanation and elaboration guideline for reporting randomised trials that use an adaptive design. BMJ, 369, m115.
6. Efron, B. (2015). Frequentist accuracy of Bayesian estimates. Journal of the Royal Statistical Society: Series B, 77(3), 617–646.
7. FDA (2026). Use of Bayesian Methodology in Clinical Trials of Drug and Biological Products: Draft Guidance for Industry. US Food and Drug Administration, January 2026.
8. Fong, E. and Holmes, C. C. (2021). Conformal Bayesian computation. Advances in Neural Information Processing Systems, 34.
9. Goligher, E. C., Tomlinson, G., Hajage, D., et al. (2018). Extracorporeal membrane oxygenation for severe acute respiratory distress syndrome and posterior probability of mortality benefit in a post hoc Bayesian analysis of a randomized clinical trial. JAMA, 320(21), 2251–2259.
10. Grünwald, P., de Heide, R. and Koolen, W. (2024). Safe testing. Journal of the Royal Statistical Society: Series B, 86(5).
11. Howard, S. R., Ramdas, A., McAuliffe, J. and Sekhon, J. (2021). Time-uniform, nonparametric, nonasymptotic confidence sequences. Annals of Statistics, 49(2), 1055–1080.
12. Johari, R., Koomen, P., Pekelis, L. and Walsh, D. (2022). Always valid inference: continuous monitoring of A/B tests. Operations Research, 70(3), 1806–1821.
13. Kohavi, R., Tang, D. and Xu, Y. (2020). Trustworthy Online Controlled Experiments: A Practical Guide to A/B Testing. Cambridge University Press.
14. Lindon, M. and Malek, A. (2022). Anytime-valid inference for multinomial count data. Advances in Neural Information Processing Systems, 35.
15. Little, R. J. (2006). Calibrated Bayes: a Bayes/frequentist roadmap. The American Statistician, 60(3), 213–223.
16. Matias, J. N., Munger, K., Le Quere, M. A. and Ebersole, C. (2021). The Upworthy Research Archive, a time series of 32,487 experiments in U.S. media. Scientific Data, 8, 195.
17. Meager, R. (2019). Understanding the average impact of microcredit expansions: a Bayesian hierarchical analysis of seven randomized experiments. American Economic Journal: Applied Economics, 11(1), 57–91.
18. Pocock, S. J. (1977). Group sequential methods in the design and analysis of clinical trials. Biometrika, 64(2), 191–199.
19. Polack, F. P., Thomas, S. J., Kitchin, N., et al. (2020). Safety and efficacy of the BNT162b2 mRNA Covid-19 vaccine. New England Journal of Medicine, 383(27), 2603–2615.
20. Ramdas, A., Grünwald, P., Vovk, V. and Shafer, G. (2023). Game-theoretic statistics and safe anytime-valid inference. Statistical Science, 38(4), 576–601.
21. REMAP-CAP Investigators (2021). Interleukin-6 receptor antagonists in critically ill patients with Covid-19. New England Journal of Medicine, 384(16), 1491–1502.
22. Robbins, H. (1970). Statistical methods related to the law of the iterated logarithm. Annals of Mathematical Statistics, 41(5), 1397–1409.
23. Robertson, D. S., Lee, K. M., López-Kolkovska, B. C. and Villar, S. S. (2023). Response-adaptive randomization in clinical trials: from myths to practical considerations. Statistical Science, 38(2), 185–208.
24. Rubin, D. B. (1984). Bayesianly justifiable and relevant frequency calculations for the applied statistician. Annals of Statistics, 12(4), 1151–1172.
25. Shafer, G. and Vovk, V. (2019). Game-Theoretic Foundations for Probability and Finance. Wiley.
26. Spiegelhalter, D. J., Abrams, K. R. and Myles, J. P. (2004). Bayesian Approaches to Clinical Trials and Health-Care Evaluation. Wiley.
27. Thompson, W. R. (1933). On the likelihood that one unknown probability exceeds another in view of the evidence of two samples. Biometrika, 25(3–4), 285–294.
28. Ville, J. (1939). Étude critique de la notion de collectif. Gauthier-Villars.
29. Wald, A. (1945). Sequential tests of statistical hypotheses. Annals of Mathematical Statistics, 16(2), 117–186.
30. Wasserman, L., Ramdas, A. and Balakrishnan, S. (2020). Universal inference. Proceedings of the National Academy of Sciences, 117(29), 16880–16890.
31. Waudby-Smith, I. and Ramdas, A. (2024). Estimating means of bounded random variables by betting. Journal of the Royal Statistical Society: Series B, 86(1), 1–27.
32. Waudby-Smith, I., Arbour, D., Sinha, R., Kennedy, E. H. and Ramdas, A. (2024). Time-uniform central limit theory and asymptotic confidence sequences. Annals of Statistics, 52(6).