Frequentists vs Bayesian power analysis

Modeling
1

Hi,

I’m trying to run a one sample bayesian power analysis using this excellent blog post.

The frequenstis version returns power of 56% while the Bayesian gives 70%.

My question is why the frequenstis power analysis differ so much from the bayesian one (especially in a simple example like this one).

As example,

I have these conditions

n_sample = 6
mu_fake = 1.87
sd_fake = 1.75

The frequentists one sample power analysis gives 56% power.

power.t.test(n = n_sample, 
                     sd =  sd_fake, 
                         delta =  mu_fake,
                         alternative = "two.sided",
                         type = "one.sample", sig.level = 0.05)

For the Bayesian I’m following the below process

# The reference model to perform one.sample t-test

# The nu prior comes from Kruschke’s [paper](https://psycnet.apa.org/fulltext/2012-18082-001.html) 
prior_t <- c(prior(normal(0, 5), class = 'Intercept'),
            prior(exponential(1), class = 'sigma'),
            prior("exponential(1.0/29)", class = 'nu'))


set.seed(917536)
bfit_reference <- brm(res ~ 1, 
            data = data.frame(res = rnorm(n_sample, mu_fake, sd_fake)), 
            family=student, 
            prior = prior_t,
            chains = 3, cores = 3, 
            silent = 2, # Mute the additional info with silent & refresh 
            refresh = 0)

# Function to get the posterior samples from the model 
post_sample <- function(x, ...)(posterior_summary(as.data.frame(x), ...))

#Simulate one new sample, update the fit and get the parameters out
sim_1s_and_fit <- function(seed, 
                   n_sample = n_sample,        
                   mu = mu_fake, # true mean
                   sd_r = sd_fake,  # repeatability error
                   model = bfit_reference
                   ) 
  {
  set.seed(seed)
    dat_f <- data.frame(res = rnorm(n_sample, mu, sd_r))
    
  updated_fit <- update(model, newdata = dat_f)
  tb_ft <- updated_fit %>% post_sample %>% as.data.frame %>% 
    rownames_to_column(var = 'param') %>% 
    filter(!grepl('lp', param)) 
  
  return(tb_ft)
}

#Apply multiple times the sim_d_and_fit() 
multiple_1s_fits_fun <- function(n_sim = 5, ...){
  tibble(seed = sample(1:1e9, n_sim)) %>% 
  mutate(b1 = map(seed, sim_1s_and_fit, ...)) %>% 
  unnest(b1) %>% return
}

(t3 <- Sys.time());
set.seed(1)
ttest_sim_ta <- multiple_1s_fits_fun(n_sim = 500, 
                                   n_sample = n_sample, mu = mu_fake, sd_r = sd_fake, 
                                   model = bfit_reference)
t4 <- Sys.time(); t4-t3

Estimate the power which in this case is around 70%

ttest_sim_ta %>% filter(param == ‘b_Intercept’) %>%
mutate(check = ifelse(Q2.5 > 0, 1, 0)) %>%
summarise(power = mean(check))

Thanks for your time.

2

I don’t have time to run your code, but I notice a few things:

  1. The Bayesian procedure you’ve outlined is analogous to checking power against a 1-tailed alternative at a 97.5% level, which is not the same as checking against a 2-tailed alternative at a 95% level. However, all else being equal I think this should result in lower apparent power as computed in your Bayesian analysis.
  2. In the Bayesian example, you’ve assumed the data are t-distributed, but in the frequentist example and in the data generation you assume the data are Gaussian. Again, I might have guessed that this would reduce the apparent power of your Bayesian analysis.
  3. You aren’t showing us your prior prior_t, which will influence your Bayesian power.
1 Like
3

Thanks for your comments, I updated my post with the priors I used.

I’m getting the 95%CI using the posterior_summary which I guess is like setting the alpha = 0.05. This returns the Q2.5 and Q97.5 so I guess this is why you say that I have the equivalant of 97.5%CI. But how could I move to the equivalent of a two tailed 95%CI?

4

The analogue of 2-tailed 95% CI would involve the fraction of fits where either the 2.5% bound is greater than zero or the 97.5% bound is less than zero.

My guess is that your priors are responsible for the bump in apparent power. What happens if you switch to a Gaussian response (equivalent to putting all the prior mass on the degrees of freedom at very large values), with a very weakly regularizing prior on the standard deviation (say, normal(0, 5))?

1 Like
5

I tried the above,

  • family=gaussian with prior for sd: normal(0, 5) the power is 0.50,
  • family=student with prior for sd: normal(0, 5) the power is 0.47

And both of them are lower than the power.t.test results.

Also I estimated the power with the below which should do the two tail 95%CI trick.

ttest_sim_ta %>% filter(param == 'b_Intercept') %>%
  mutate(check = ifelse(Q2.5 > 0 | Q97.5 < 0, 1, 0)) %>% 
  summarise(power =  mean(check))

The prior seems to have a nice regularizing effect in this simple power analysis example.
Thanks for your help, I think this is fascinating.

6

This is a somewhat recurrent and possibly controversial discussion, but technically the so called “Frequentist” approaches are equivalent to Bayesian inference up to the specification of priors, and therefore the former are a particular case of the latter, strictly speaking.

In practice, there end up being several differences between “Frequentists” and “Bayesians” beyond priors, but those are differences of practice, not of theoretical underpinnings. For your particular case @jsocolar points out clearly the equivalence and how to replicate it in a Bayesian framework, and in general if you understand the formal “Frequentists” statistical description it is straightforward to formulate a “Bayesian” version.

2 Likes