A little philosophy

Any post describing Bayesian models, p-values, and Type I errors – even one focused on simulation – would be remiss without at least some discussion of underlying philosophical differences between the two approaches. To start, this statement succinctly describes a key conflict:

The Bayesian approach is based on determining the probability of a hypothesis with a model using an “a priori” probability that is then updated based on data. On the contrary, [in] classical hypothesis testing … the error-I type is the probability of wrongly refusing the null hypothesis when it is true. Thus, [frequentist logic] is completely different from Bayesian logic (since probability is referred to making a mistake, not to the hypothesis itself).

So, the Bayesian approach is concerned with estimating \(P(hypothesis \ | \ data)\), whereas the frequentist p-value is an estimate of \(P(data \ | \ null \ hypothesis)\), quite different animals and it isn’t obvious they can be reconciled. Indeed, as Frank Harrell points out, there is a logical inconsistency with assuming that a particular (null) hypothesis is true and then trying to draw conclusions about that very same hypothesis:

The probability of a treatment not being effective is the probability of “regulator’s regret.’’ One must be very clear on what is conditioned upon (assumed) in computing this probability. Does one condition on the true effectiveness or does one condition on the available data? Type I error conditions on the treatment having no effect … Can one quantify evidence for making a wrong decision if one assumes that all conclusions of non-zero effect are wrong up front because \(H_0\) was assumed to be true? Aren’t useful error probabilities the ones that are not based on assumptions about what we are assessing but rather just on the data available to us? [emphasis added]

In addition to Harrell’s post, there is a paper by Greenland et al that describes how statistical tests, confidence intervals, and statistical power can easily be misinterpreted and abused. However, these points of view have not been fully adopted by the scientific community more broadly, and Type I error rates and statistical power are still the predominant way of evaluating study designs.

This paper by Ryan et al acknowledges these realities while asking whether “we need to adjust for interim analyses in a Bayesian adaptive trial design?” The conclusion is “yes” if the goal is indeed to satisfy operating characteristics defined within the frequentist paradigm – and I will get to that in a second – but they conclude, maybe a little wistfully, that

if we avoid this dichotomisation [reject or fail to reject] and simply report the posterior probability of benefit, then we could potentially avoid having to specify the type I error of a Bayesian design.

But until more reviewers and funders accept this, investigators who opt for a Bayesian approach will likely still be required to square this circle.

Simulating interim looks

To assess the Type I error, I’ve used a relatively simple data generating process and estimation model. There is a continuous outcome \(Y\) that is normally distributed, and the mean is entirely a function of the treatment arm assignment \(Z\). When \(Z_i = 0\), the subject is in the control arm and the mean is \(\alpha\); when \(Z_i = 1\), the subject is in the treatment arm, and the mean is \(\alpha + \beta\). The standard deviation for both groups is assumed to be the same, \(\sigma_s\):

\[ Y_i \sim N(\alpha + \beta Z_i, \sigma_s) \]

To fit the model, we have to assume prior distributions for the parameters \(\alpha\), \(\beta\), and \(\sigma_s\); here is the assumption for \(\beta\), which is the treatment effect, the parameter of primary interest:

\[ \beta \sim t_\text{student}(df = 3, \mu = 0, \sigma = 10) \]

In the model specification below, I have allowed the prior standard deviation of \(\beta\) to be specified as an argument, so that we can explore the impact of different assumptions on the operating characteristics. Here is the Stan code to implement the model:

data {
  int<lower=0> N;
  int<lower=0,upper=1> rx[N];
  vector[N] y;
  real p_mu;
  real p_sigma;
}

parameters {
  real alpha;
  real beta;
  real<lower=0> sigma;
}

transformed parameters {
  real yhat[N];
  
  for (i in 1:N) 
    yhat[i] = alpha + beta * rx[i];
}

model {
  alpha ~ student_t(3, 0, 10);
  beta ~ student_t(3, p_mu, p_sigma);
  
  y ~ normal(yhat, sigma);
}

And here is a call to compile the code using package cmdstanr:

library(simstudy)
library(data.table)
library(cmdstanr)
library(posterior)

mod <- cmdstan_model("code/multiple.stan")

I’ve written a function bayes_fit that (1) estimates the model, (2) collects the samples of \(\beta\) from the posterior distribution, and (3) returns the probability of “success.” Here, success is determined one of two ways. In the first, we evaluate \(P(\beta > 0)\) (using the posterior distribution) and declare success if this probability exceeds 95%; this is the probability that the intervention is successful. The second decision rule uses the first criterion and adds an additional one to ensure that the effect size is clinically meaningful by requiring \(P(\beta > M) > 50\%\). \(M\) is some meaningful threshold that has been agreed upon prior to conducting the study. Under the two-part decision rule, both requirements ( \(P(\beta > 0) > 95\%\) and \(P(\beta > M) > 50\%\)) need to be satisfied in order to declare success.

bayes_fit <- function(dx, p_sigma, m_effect, decision_rule, x) {
  
  # 1: estimate model
  
  data_list <- list(N = nrow(dx), y = dx$y, rx = dx$rx, p_mu = 0, p_sigma = p_sigma)

  fit <- mod$sample(
    data = data_list,
    refresh = 0,
    chains = 4L,
    parallel_chains = 4L,
    iter_warmup = 500,
    iter_sampling = 2500,
    step_size = 0.1,
    show_messages = FALSE
  )
  
  # 2: collect sample of betas from posterior
  
  df <- data.frame(as_draws_rvars(fit$draws(variables = "beta")))
  
  # 3: evaluate success based on desired decision rule
  
  if (decision_rule == 1) {
    return((mean(df$beta > 0) > 0.95))
  } else { # decision_rule == 2
    return( ((mean(df$beta > 0) > 0.95) & (mean(df$beta > m_effect ) > 0.5)) )  
  }
}

Function freq_fit fits a linear model and returns a p-value for the estimate of \(\beta\), to be used as a basis for comparison with the Bayesian models:

freq_fit <- function(dx) {
  
  lmfit <- lm(y ~ rx, data = dx)
  coef(summary(lmfit))["rx", "Pr(>|t|)"]
  
}

The data generation process is simple and assumes that \(\alpha = 0\), and more importantly, that the treatment effect \(\beta=0\) (and that \(\sigma_s = 1\))

def <- defData(varname = "rx", formula = "1;1", dist = "trtAssign")
def <- defData(def, varname = "y", formula = 0, variance = 1, dist = "normal")
  
set.seed(1918721)
dd <- genData(1000, def)

dd
##         id rx      y
##    1:    1  1 -0.593
##    2:    2  1  0.648
##    3:    3  0 -0.819
##    4:    4  0 -0.386
##    5:    5  0 -0.304
##   ---               
##  996:  996  1  0.343
##  997:  997  0 -0.948
##  998:  998  1 -0.497
##  999:  999  0 -1.301
## 1000: 1000  0  0.282

And here we are at the key point where we simulate the interim looks. I’ve taken a crude, probably inefficient approach that has the advantage of being extremely easy to code. Using the R function sapply, we can sequentially estimate models using incrementally larger segments of a particular data set. In the first round, the model is estimated for the first 100 observations, in the second round the model is estimated for the first 200, and so on. For each round, we get an indicator for whether the trial would have been declared successful based on whatever criteria were being used. If any of the interim looks results in a success, then that simulated study would be deemed (inappropriately) successful. This approach is inherently inefficient, because we are conducting all interim looks regardless of the outcome of earlier looks; in effect, we are fitting too many models. But in the case where only a relatively small number of studies will be “successful,” this is a small price to pay for ease of coding.

bayes_ci <- sapply(seq(100, 1000, by = 100), 
    function(x) bayes_fit(dd[1:x], p_sigma = 10, m_effect = 0.2, x))

The function returns a vector of success indicators, one for each interim look.

bayes_ci
##  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

Here is a visualization representing the underlying posterior densities at each interim look. It is clear that less than 95% of each density (shaded in darker red) falls to the right of \(0\) at each stage of evaluation, and based on the single decision rule, the study would not have been declared successful at any interim look:

In practice, all we need to do is check whether any of the looks has resulted in a “success”. In this case, there were no successes:

any(bayes_ci)
## [1] FALSE

Operating characteristics: Type I error rates

##    p_sigma m_effect  freq  bayes
## 1:       1        0 0.118 0.0937
## 2:       5        0 0.099 0.0993
## 3:      10        0 0.101 0.0917

##    p_sigma m_effect  freq bayes
## 1:       1        0 0.193 0.157
## 2:       5        0 0.197 0.166
## 3:      10        0 0.189 0.169

Addendum

Code used to generate p-value estimates for a variety of scenarios was written to be executed on a high performance computing cluster, which is a highly parallel computing environment. I am using package rslurm to manage the parallelization of the process (there is an alternative package called slurmR which I have used in the past). The replications of each of the scenarios (represented by a set of simulation parameters) are stored in data frame scenarios, and the function slurm_apply cycles through the data frame, with many replications executed in parallel. The function s_replication calls the two functions bayes_fit and freq_fit described above.

library(simstudy)
library(data.table)
library(parallel)
library(cmdstanr)
library(posterior)
library(rslurm)

### Function to generate data and estimate parameters

s_replicate <- function(iter, p_sigma, decision_rule, m_effect, seq) {
  
  set_cmdstan_path(path = "/gpfs/.../cmdstan/2.25.0")
  
  def <- defData(varname = "rx", formula = "1;1", dist = "trtAssign")
  def <- defData(def, varname = "y", formula = 0, variance = 1, dist = "normal")
  
  dd <- genData(end, def)
  
  freq_ps <- sapply(seq(start, end, by = by), function(x) freq_fit(dd[1:x]))
  freq_effect <- any(freq_ps < 0.05)
  
  bayes_ci <- sapply(seq(start, end, by = by), 
    function(x) bayes_fit(dd[1:x], p_sigma, m_effect, decision_rule, x))
  bayes_effect <- any(bayes_ci)
  
  return(data.table(seq, iter, p_sigma, m_effect, decision_rule, 
    freq_effect, bayes_effect))  
}

### Set simulation parameters

scenario_dt <- function(...) {
  argdt <- data.table(expand.grid(...))
  argdt[, seq := .I]
  argdt[]
}

iter <- c(1:1000)
p_sigma <- c(1, 5, 10)
decision_rule = 2
m_effect <- c(0.2, 0.3, 0.4) # if decision_rule = 2
# decision_rule = 1
# m_effect <- 0

start <- 100L
end <- 1000L
by <- 100L

scenarios <- scenario_dt(
  iter = iter, 
  p_sigma = p_sigma, 
  decision_rule = decision_rule,
  m_effect = m_effect
)

### Compile stan code

set_cmdstan_path(path = "/gpfs/.../cmdstan/2.25.0")
mod <- cmdstan_model("multiple.stan")

### Set rslurm arguments

sopts <- list(time = '12:00:00', partition = "cpu_short", `mem-per-cpu` = "5G")
sobjs <- c("freq_fit", "bayes_fit", "mod", "start", "end", "by")

### Replicate over iterations

sjob <- slurm_apply(
  f = s_replicate, # the function
  params = scenarios, # a data frame
  jobname = 'mult_i',
  nodes = 50, 
  slurm_options = sopts,
  global_objects = sobjs,
  submit = TRUE
)

### Collect the results and save them

res <- get_slurm_out(sjob, outtype = 'table', wait = TRUE)
save(res, file = "/gpfs/.../mult.rda")