Filtering by Significance Mis-estimates Small Effect Sizes

This is a demonstration of what has been pointed out by Gelman and colleagues that when dealing with truly small effects, screening results by statistical significance will leade to large Type-M (magnitude) and Type-S (sign) errors.

I’ve actually blogged about this before, but I’m experimenting with using tidyverse functions to do the actual simulation.

library(tidyverse)

Loading tidyverse: ggplot2
Loading tidyverse: tibble
Loading tidyverse: tidyr
Loading tidyverse: readr
Loading tidyverse: purrr
Loading tidyverse: dplyr
Conflicts with tidy packages -------------------------------------------------
filter(): dplyr, stats
lag():    dplyr, stats

library(broom)
library(scales)

Attaching package: ‘scales’

The following object is masked from ‘package:purrr’:

    discard

The following object is masked from ‘package:readr’:

    col_factor

The Effect Size

Type-M errors are most likely to happen when the true effect is small. I’ll be looking at the difference between normal distributions with \(\mu=1\) and \(\mu=1.25\) and \(\sigma = 1\), which would constitute a Cohen’s d of 0.25.

group_a_dens <- function(x){
  dnorm(x, mean = 1, sd = 1)
}

group_b_dens <- function(x){
  dnorm(x, mean = 1.25, sd = 1)
}

dummy_df <- data.frame(x = 1, y = 1)
ggplot(dummy_df, aes(x, y))+
  stat_function(fun = group_a_dens, 
                geom = "line",
                aes(color = 'group_a'))+
  stat_function(fun = group_b_dens, 
                geom = "line",
                aes(color = 'group_b'))+  
  xlim(-2,4.25)+
  scale_color_brewer("group", palette = "Dark2")+
  theme_minimal()

The Simulation

I’ll simulate samples of various sizes (N = 10, 20, 50, 100, 500, 1,000) from these two groups, 1,000 times for each sample size. I’ll perform a t-test for each simulated sample.

make_samp <- function(x, N){
  data_frame(iter = x,
             group_a = rnorm(N, mean = 1, sd = 1),
             group_b = rnorm(N, mean = 1.25, sd = 1))
}

do_t_test <- function(df){
  t.test(df$group_a, df$group_b)
}

This sets up the simulation. It creates one row for each simulation.

sim_df <- expand.grid(iter = 1:1000,
                      size = c(10, 20, 50, 100, 500, 1000))%>%
          tbl_df()

The real trick here is using map2() to run it.

sims <- sim_df %>%
            mutate(sims = map2(iter, size, make_samp),
                   t_test = map(sims, do_t_test),
                   params = map(t_test, tidy))%>%
            unnest(params)

sim_agg <- sims %>%
            group_by(size, 
                     sig = p.value <= 0.05,
                     sign = sign(estimate))%>%
            summarise(effect = mean(estimate),
                      lo = quantile(estimate, 0.025),
                      hi = quantile(estimate, 0.975),
                      N = n())%>%
            filter(sig)

After filtering out the non-significant tests, I’ll be plotting the estimated effect/-0.25 (the true effect size) to show how many times larger the it than the true effect.

sim_agg %>%
  ggplot(aes(factor(size), effect/-0.25))+
    geom_hline(yintercept = 1)+
    geom_point(aes(size = N/1000))+
    geom_segment(aes(y = lo/-0.25, 
                     yend = hi/-0.25,
                     xend = factor(size)))+
    theme_minimal()+
    ggtitle("p <= 0.05")

The smaller sample sizes have a high Type 1 error rate because they are underpowered. However, the average estimated effect size for the significant tests are highly over-estimated. This over-estimation is not ameliorated by using a different p-value threshold (that just boosts the rate of Type 1 errors).

sim_agg_01 <- sims %>%
            group_by(size, 
                     sig = p.value <= 0.01,
                     sign = sign(estimate))%>%
            summarise(effect = mean(estimate),
                      lo = quantile(estimate, 0.025),
                      hi = quantile(estimate, 0.975),
                      N = n())%>%
            filter(sig)

sim_agg_01 %>%
  ggplot(aes(factor(size), effect/-0.25))+
    geom_hline(yintercept = 1)+
    geom_point(aes(size = N/1000))+
    geom_segment(aes(y = lo/-0.25, 
                     yend = hi/-0.25,
                     xend = factor(size)))+
    theme_minimal()+
    ggtitle("p <= 0.01")