Bayesian_power_analysis

Bayesian_power_analysis.R

# |- table.express ----
library(tidyverse)
library(data.table)
library(brms)
library(broom)

# plot theme
theme_set(theme_dark() +
            theme(legend.position = "none",
                  panel.grid = element_blank(),
                  panel.background = element_rect(fill = "grey20"),
                  plot.background = element_rect(fill = "grey20"),
                  text = element_text(color = "white"),
                  axis.text = element_text(color = "white")))

# data
d <- data.table(yc = rnorm(100, 0, 1),
                yt = rnorm(100, .5, 1),
                id = 1:100) %>% 
  melt(id.vars = "id", measure.vars = c("yc","yt"), variable.name = "group")

# model
fit <- brm(data = d,
           family = gaussian,
           value ~ 0 + intercept + group,
           prior = c(prior(normal(0, 10), class = b),
                     prior(student_t(3, 1, 10), class = sigma)),
           seed = 1,
           chains = 1)

library(table.express)

# function to simulate new data for model
sim_d <- function(seed, n) {
  set.seed(seed)
  
  mu_t <- .5
  mu_c <- 0
  
  data.table(group = rep(c("control", "treatment"), n/2)) %>% 
    mutate(value = ifelse(group == "control",
                          rnorm(n, mean = mu_c, sd = 1), 
                          rnorm(n, mean = mu_t, sd = 1)))
}

n_sim <- 100
n <- 80

# single simulation
sims <- data.table(seed = 1:n_sim) %>% 
  mutate(d = map(seed, sim_d, n = n)) %>% 
  mutate(fit = map2(d, seed, ~update(fit, newdata = .x, seed = .y)))

# tidy and extract model estimates
sims <- sims %>% 
  mutate(effects = map(fit, ~tidy(.x, prob = .95))) %>% 
  group_by(seed) %>% 
  select(unlist(effects, recursive = FALSE)) %>% 
  filter(term == "b_grouptreatment") %>% 
  mutate(powered = case_when(lower > 0 ~ 1, TRUE ~ 0))

sims %>% 
  group_by(powered) %>% 
  select(.N)

# multiple simulations
cur <- tibble(seed = 1:n_sim) %>% 
  crossing(sample = seq(40, 140, by = 20)) %>% 
  as.data.table() %>% 
  mutate(d = map2(seed, sample, ~sim_d(seed = .x, n = .y))) %>% 
  mutate(fit = map2(d, seed, ~update(fit, newdata = .x, seed = .y))) %>% 
  mutate(effects = map(fit, ~tidy(.x, prob = .95))) %>% 
  group_by(seed, sample) %>% 
  select(unlist(effects, recursive = FALSE)) %>% 
  filter(term == "b_grouptreatment") %>% 
  mutate(powered = case_when(lower > 0 ~ 1, TRUE ~ 0))

# plots
cur %>% 
  group_by(sample) %>% 
  select(mean(powered)) %>% 
  ggplot(aes(sample, V1)) +
  geom_point(color = "grey80") +
  geom_line(color = "grey80") +
  geom_hline(yintercept = .8, color = "grey80", linetype = "dashed") +
  labs(x = "Sample Size", 
       y = "Power")

ggplot(cur, aes(seed, estimate, ymin = lower, ymax = upper)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey70") +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "grey70") +
  geom_pointrange(size = .1, aes(color = factor(powered))) +
  labs(y = "Estimate (Credibility Interval)",
       x = "Index") +
  scale_color_manual(values = c("dodgerblue1", "#F1C40F")) +
  facet_wrap(~sample)