Skip to content

Instantly share code, notes, and snippets.

@dgrtwo
Last active May 11, 2019 20:01
Show Gist options
  • Save dgrtwo/1337b9d5c80adde1bc7a72048d9e4cc6 to your computer and use it in GitHub Desktop.
Save dgrtwo/1337b9d5c80adde1bc7a72048d9e4cc6 to your computer and use it in GitHub Desktop.
# 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.")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment