underpowered-simulation.R

# Code behind this tweet:

# https://twitter.com/drob/status/1126988304090574848

library(tidyverse)
library(broom)

t_tests <- crossing(pi0 = .75,
                    effect_size = .25,
                    trial = 1:10000,
                    sample_size = c(10, 30, 100, 300)) %>%
  mutate(h1 = runif(n()) > pi0) %>%
  unnest(map(sample_size, seq_len)) %>%
  mutate(x = rnorm(n(), h1 * effect_size, 1)) %>%
  group_by(pi0, effect_size, sample_size, trial, h1) %>%
  summarize(test = list(t.test(x))) %>%
  unnest(map(test, tidy)) %>%
  ungroup()

t_tests %>%
  ggplot(aes(p.value, fill = ifelse(h1, "alternative", "null"))) +
  geom_histogram(binwidth = .05, boundary = 0) +
  facet_wrap(~ sample_size, scales = "free_y") +
  labs(fill = "H1",
       x = "P-value",
       title = "How p-value distribution changes with sample size",
       subtitle = "One-sample t-test, pi0 = .75, alternative mean = .25, m = 10000")

# And followup:

library(qvalue)

qvals <- t_tests %>%
  group_by(sample_size) %>%
  mutate(qvalue = qvalue(p.value)$qvalues) %>%
  summarize(pi1_hat = 1 - qvalue(p.value)$pi0,
            discoveries = mean(qvalue < .05))

qvals %>%
  gather(metric, value, -sample_size) %>%
  mutate(metric = ifelse(metric == "discoveries",
                         "% significant at .05 FDR",
                         "Estimated π_1")) %>%
  ggplot(aes(sample_size, value, color = metric)) +
  geom_line() +
  geom_hline(yintercept = .25, lty = 2) +
  scale_x_log10() +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(x = "Sample size",
       y = "Percentage of hypotheses",
       color = "",
       title = "Estimating how underpowered a test is",
       subtitle = "One-sample t-test, alternative mean = .25, m = 10000. True π_1 shown as dashed line.")