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")

LS0tCnRpdGxlOiAiRmlsdGVyaW5nIGJ5IFNpZ25pZmljYW5jZSBNaXMtZXN0aW1hdGVzIFNtYWxsIEVmZmVjdCBTaXplcyIKYXV0aG9yOiAiSm9zZWYgRnJ1ZWh3YWxkIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpUaGlzIGlzIGEgZGVtb25zdHJhdGlvbiBvZiB3aGF0IGhhcyBiZWVuIHBvaW50ZWQgb3V0IGJ5IFtHZWxtYW4gYW5kIGNvbGxlYWd1ZXNdKGh0dHA6Ly9hbmRyZXdnZWxtYW4uY29tLzIwMTAvMTIvMTMvdGhlX3RydXRoX3dlYXJzLykgdGhhdCB3aGVuIGRlYWxpbmcgd2l0aCB0cnVseSBzbWFsbCBlZmZlY3RzLCBzY3JlZW5pbmcgcmVzdWx0cyBieSBzdGF0aXN0aWNhbCBzaWduaWZpY2FuY2Ugd2lsbCBsZWFkZSB0byBsYXJnZSBUeXBlLU0gKG1hZ25pdHVkZSkgYW5kIFR5cGUtUyAoc2lnbikgZXJyb3JzLgoKSSd2ZSBhY3R1YWxseSBibG9nZ2VkIGFib3V0IHRoaXMgYmVmb3JlLCBidXQgSSdtIGV4cGVyaW1lbnRpbmcgd2l0aCB1c2luZyBgdGlkeXZlcnNlYCBmdW5jdGlvbnMgdG8gZG8gdGhlIGFjdHVhbCBzaW11bGF0aW9uLgoKYGBge3J9CmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KGJyb29tKQpsaWJyYXJ5KHNjYWxlcykKYGBgCgojIFRoZSBFZmZlY3QgU2l6ZQoKVHlwZS1NIGVycm9ycyBhcmUgbW9zdCBsaWtlbHkgdG8gaGFwcGVuIHdoZW4gdGhlIHRydWUgZWZmZWN0IGlzIHNtYWxsLiBJJ2xsIGJlIGxvb2tpbmcgYXQgdGhlIGRpZmZlcmVuY2UgYmV0d2VlbiBub3JtYWwgZGlzdHJpYnV0aW9ucyB3aXRoICRcbXU9MSQgYW5kICRcbXU9MS4yNSQgYW5kICRcc2lnbWEgPSAxJCwgd2hpY2ggd291bGQgY29uc3RpdHV0ZSBhIENvaGVuJ3MgZCBvZiAwLjI1LiAKCgpgYGB7ciB0aWR5ID0gRn0KZ3JvdXBfYV9kZW5zIDwtIGZ1bmN0aW9uKHgpewogIGRub3JtKHgsIG1lYW4gPSAxLCBzZCA9IDEpCn0KYGBgCgpgYGB7cn0KZ3JvdXBfYl9kZW5zIDwtIGZ1bmN0aW9uKHgpewogIGRub3JtKHgsIG1lYW4gPSAxLjI1LCBzZCA9IDEpCn0KYGBgCgpgYGB7cn0KZHVtbXlfZGYgPC0gZGF0YS5mcmFtZSh4ID0gMSwgeSA9IDEpCgpnZ3Bsb3QoZHVtbXlfZGYsIGFlcyh4LCB5KSkrCiAgc3RhdF9mdW5jdGlvbihmdW4gPSBncm91cF9hX2RlbnMsIAogICAgICAgICAgICAgICAgZ2VvbSA9ICJsaW5lIiwKICAgICAgICAgICAgICAgIGFlcyhjb2xvciA9ICdncm91cF9hJykpKwogIHN0YXRfZnVuY3Rpb24oZnVuID0gZ3JvdXBfYl9kZW5zLCAKICAgICAgICAgICAgICAgIGdlb20gPSAibGluZSIsCiAgICAgICAgICAgICAgICBhZXMoY29sb3IgPSAnZ3JvdXBfYicpKSsgIAogIHhsaW0oLTIsNC4yNSkrCiAgc2NhbGVfY29sb3JfYnJld2VyKCJncm91cCIsIHBhbGV0dGUgPSAiRGFyazIiKSsKICB0aGVtZV9taW5pbWFsKCkKYGBgCgoKCgoKIyBUaGUgU2ltdWxhdGlvbgoKSSdsbCBzaW11bGF0ZSBzYW1wbGVzIG9mIHZhcmlvdXMgc2l6ZXMgKE4gPSAxMCwgMjAsIDUwLCAxMDAsIDUwMCwgMSwwMDApIGZyb20gdGhlc2UgdHdvIGdyb3VwcywgMSwwMDAgdGltZXMgZm9yIGVhY2ggc2FtcGxlIHNpemUuIEknbGwgcGVyZm9ybSBhIHQtdGVzdCBmb3IgZWFjaCBzaW11bGF0ZWQgc2FtcGxlLgoKYGBge3J9Cm1ha2Vfc2FtcCA8LSBmdW5jdGlvbih4LCBOKXsKICBkYXRhX2ZyYW1lKGl0ZXIgPSB4LAogICAgICAgICAgICAgZ3JvdXBfYSA9IHJub3JtKE4sIG1lYW4gPSAxLCBzZCA9IDEpLAogICAgICAgICAgICAgZ3JvdXBfYiA9IHJub3JtKE4sIG1lYW4gPSAxLjI1LCBzZCA9IDEpKQp9CmBgYAoKYGBge3J9Cgpkb190X3Rlc3QgPC0gZnVuY3Rpb24oZGYpewogIHQudGVzdChkZiRncm91cF9hLCBkZiRncm91cF9iKQp9CmBgYAoKClRoaXMgc2V0cyB1cCB0aGUgc2ltdWxhdGlvbi4gSXQgY3JlYXRlcyBvbmUgcm93IGZvciBlYWNoIHNpbXVsYXRpb24uIAoKYGBge3J9CnNpbV9kZiA8LSBleHBhbmQuZ3JpZChpdGVyID0gMToxMDAwLAogICAgICAgICAgICAgICAgICAgICAgc2l6ZSA9IGMoMTAsIDIwLCA1MCwgMTAwLCA1MDAsIDEwMDApKSU+JQogICAgICAgICAgdGJsX2RmKCkKCmBgYAoKVGhlIHJlYWwgdHJpY2sgaGVyZSBpcyB1c2luZyBgbWFwMigpYCB0byBydW4gaXQuCgpgYGB7cn0Kc2ltcyA8LSBzaW1fZGYgJT4lCiAgICAgICAgICAgIG11dGF0ZShzaW1zID0gbWFwMihpdGVyLCBzaXplLCBtYWtlX3NhbXApLAogICAgICAgICAgICAgICAgICAgdF90ZXN0ID0gbWFwKHNpbXMsIGRvX3RfdGVzdCksCiAgICAgICAgICAgICAgICAgICBwYXJhbXMgPSBtYXAodF90ZXN0LCB0aWR5KSklPiUKICAgICAgICAgICAgdW5uZXN0KHBhcmFtcykKYGBgCgpgYGB7cn0Kc2ltX2FnZyA8LSBzaW1zICU+JQogICAgICAgICAgICBncm91cF9ieShzaXplLCAKICAgICAgICAgICAgICAgICAgICAgc2lnID0gcC52YWx1ZSA8PSAwLjA1LAogICAgICAgICAgICAgICAgICAgICBzaWduID0gc2lnbihlc3RpbWF0ZSkpJT4lCiAgICAgICAgICAgIHN1bW1hcmlzZShlZmZlY3QgPSBtZWFuKGVzdGltYXRlKSwKICAgICAgICAgICAgICAgICAgICAgIGxvID0gcXVhbnRpbGUoZXN0aW1hdGUsIDAuMDI1KSwKICAgICAgICAgICAgICAgICAgICAgIGhpID0gcXVhbnRpbGUoZXN0aW1hdGUsIDAuOTc1KSwKICAgICAgICAgICAgICAgICAgICAgIE4gPSBuKCkpJT4lCiAgICAgICAgICAgIGZpbHRlcihzaWcpCmBgYAoKQWZ0ZXIgZmlsdGVyaW5nIG91dCB0aGUgbm9uLXNpZ25pZmljYW50IHRlc3RzLCBJJ2xsIGJlIHBsb3R0aW5nIHRoZSBlc3RpbWF0ZWQgZWZmZWN0Ly0wLjI1ICh0aGUgdHJ1ZSBlZmZlY3Qgc2l6ZSkgdG8gc2hvdyBob3cgbWFueSB0aW1lcyBsYXJnZXIgdGhlIGl0IHRoYW4gdGhlIHRydWUgZWZmZWN0LgoKCgpgYGB7cn0Kc2ltX2FnZyAlPiUKICBnZ3Bsb3QoYWVzKGZhY3RvcihzaXplKSwgZWZmZWN0Ly0wLjI1KSkrCiAgICBnZW9tX2hsaW5lKHlpbnRlcmNlcHQgPSAxKSsKICAgIGdlb21fcG9pbnQoYWVzKHNpemUgPSBOLzEwMDApKSsKICAgIGdlb21fc2VnbWVudChhZXMoeSA9IGxvLy0wLjI1LCAKICAgICAgICAgICAgICAgICAgICAgeWVuZCA9IGhpLy0wLjI1LAogICAgICAgICAgICAgICAgICAgICB4ZW5kID0gZmFjdG9yKHNpemUpKSkrCiAgICB0aGVtZV9taW5pbWFsKCkrCiAgICBnZ3RpdGxlKCJwIDw9IDAuMDUiKQpgYGAKCgpUaGUgc21hbGxlciBzYW1wbGUgc2l6ZXMgaGF2ZSBhIGhpZ2ggVHlwZSAxIGVycm9yIHJhdGUgYmVjYXVzZSB0aGV5IGFyZSB1bmRlcnBvd2VyZWQuIEhvd2V2ZXIsIHRoZSBhdmVyYWdlIGVzdGltYXRlZCBlZmZlY3Qgc2l6ZSBmb3IgdGhlIHNpZ25pZmljYW50IHRlc3RzIGFyZSBoaWdobHkgb3Zlci1lc3RpbWF0ZWQuIFRoaXMgb3Zlci1lc3RpbWF0aW9uIGlzIG5vdCBhbWVsaW9yYXRlZCBieSB1c2luZyBhIGRpZmZlcmVudCBwLXZhbHVlIHRocmVzaG9sZCAodGhhdCBqdXN0IGJvb3N0cyB0aGUgcmF0ZSBvZiBUeXBlIDEgZXJyb3JzKS4KCmBgYHtyfQpzaW1fYWdnXzAxIDwtIHNpbXMgJT4lCiAgICAgICAgICAgIGdyb3VwX2J5KHNpemUsIAogICAgICAgICAgICAgICAgICAgICBzaWcgPSBwLnZhbHVlIDw9IDAuMDEsCiAgICAgICAgICAgICAgICAgICAgIHNpZ24gPSBzaWduKGVzdGltYXRlKSklPiUKICAgICAgICAgICAgc3VtbWFyaXNlKGVmZmVjdCA9IG1lYW4oZXN0aW1hdGUpLAogICAgICAgICAgICAgICAgICAgICAgbG8gPSBxdWFudGlsZShlc3RpbWF0ZSwgMC4wMjUpLAogICAgICAgICAgICAgICAgICAgICAgaGkgPSBxdWFudGlsZShlc3RpbWF0ZSwgMC45NzUpLAogICAgICAgICAgICAgICAgICAgICAgTiA9IG4oKSklPiUKICAgICAgICAgICAgZmlsdGVyKHNpZykKYGBgCgpgYGB7cn0Kc2ltX2FnZ18wMSAlPiUKICBnZ3Bsb3QoYWVzKGZhY3RvcihzaXplKSwgZWZmZWN0Ly0wLjI1KSkrCiAgICBnZW9tX2hsaW5lKHlpbnRlcmNlcHQgPSAxKSsKICAgIGdlb21fcG9pbnQoYWVzKHNpemUgPSBOLzEwMDApKSsKICAgIGdlb21fc2VnbWVudChhZXMoeSA9IGxvLy0wLjI1LCAKICAgICAgICAgICAgICAgICAgICAgeWVuZCA9IGhpLy0wLjI1LAogICAgICAgICAgICAgICAgICAgICB4ZW5kID0gZmFjdG9yKHNpemUpKSkrCiAgICB0aGVtZV9taW5pbWFsKCkrCiAgICBnZ3RpdGxlKCJwIDw9IDAuMDEiKQpgYGAK