This is the script which is mentioned in the article HERE.

library(ggplot2) 

sim = 1000 

x1 = y1 = NULL 

for (n in c(16,32,64))  

{ 

    p.values_all = p.values_sex = p.values_int = p.values_m = p.values_f =  

    p.values_rep = p.values_neg = p.values_pos = NULL 

    for (iteration in 1:sim) 

  { 

    p.value_now = NULL 

    all1 = rnorm(n) 

    all2 = rnorm(n) 

    m1 = all1[1:(n/2)] 

    m2 = all2[1:(n/2)] 

    f1 = all1[((n/2)+1):n] 

    f2 = all2[((n/2)+1):n] 

    data1 <- c(all1,all2) 

    sex1 <- c(rep(c('M', 'F'), n)) 

    group1 <- c(rep('A', n),rep('B', n)) 

    sabv <- data.frame(group1,sex1,data1) 

    resu <- aov(data1 ~ group1 + sex1+ group1*sex1, data=sabv) 

    p.value_now <- t.test(all1,all2)$p.value 

    if (p.value_now > 0.05){ 

      p.values_neg <- p.value_now 

      p.value_now <- t.test(m1,m2)$p.value 

      if (p.value_now > 0.05){ 

        p.value_now <- t.test(f1,f2)$p.value 

        if (p.value_now > 0.05){ 

          p.value_now <- summary(resu)[[1]][[2,"Pr(>F)"]] 

          if (p.value_now > 0.05){ 

            p.value_now <- summary(resu)[[1]][[3,"Pr(>F)"]] 

            if (p.value_now > 0.05){ 

              p.values_rep = c(p.values_rep, p.values_neg) 

              next 

            }   

            p.values_int = c(p.values_int, summary(resu)[[1]][[3,"Pr(>F)"]]) 

            next 

          }    

          p.values_sex = c(p.values_sex, summary(resu)[[1]][[2,"Pr(>F)"]]) 

          next 

        }   

        p.values_f = c(p.values_f, t.test(f1,f2)$p.value) 

        next 

      } 

      p.values_m = c(p.values_m, t.test(m1,m2)$p.value) 

      next 

    } 

    p.values_pos <- c(p.values_pos, p.value_now) 

    next 

  } 

  y1 <- c(y1, c(p.values_m, p.values_f, p.values_rep, p.values_int)) 

  yn <- length(c(p.values_m, p.values_f, p.values_rep, p.values_int)) 

  x1 <- c(x1, c(rep(n,yn))) 

} 

fickle <- data.frame(x1,y1) 

fickle_plot <- ggplot(fickle, aes(x=x1, y=y1)) 

fickle_plot + xlab("Sample size") + ylab("P values") +  

  geom_jitter(shape=16, size=2, width = 0.06) + 

  scale_y_log10(breaks=c(0.00001,0.0001,0.001,0.01,0.05,0.1,1), labels = scales::comma) + 

  scale_x_log10(breaks=c(16,32,64)) + 

  theme(panel.grid.major=element_blank()) + 

  geom_hline(yintercept=0.05, linetype="dashed", color = "red") + 

  theme(axis.title.y = element_text(margin = margin(r = 15))) + 

  theme(axis.title.x = element_text(margin = margin(t = 15))) + 

  theme(plot.margin = margin(1, 1, 1, 1, "cm"))