sim = function(n1       = 10,
               n2       = 10,
               delta    = 0,
               sd       = 1,
               alpha    = .05,
               chk_norm = TRUE,
               norm_a   = 0.05,
               rem_out  = TRUE,
               out_a    = 0.1){

  ad2 = alpha/2
  # Gen prng data from normal dists
  a = rnorm(n1, 0, sd)
  b = rnorm(n2, delta, sd)


  # If checkign for outliers filter any out if outside the specified percentile
  if(rem_out){
    a_out = abs(.5 - ecdf(a)(a)) < ad2
    b_out = abs(.5 - ecdf(b)(b)) < ad2
    if(mean(a_out) >= out_a | mean(b_out) >= out_a){
      a = a[!a_out]
      b = b[!b_out]
    }
  }

  # Do normality test for both groups
  norm_pA = shapiro.test(a)$p.value
  norm_pB = shapiro.test(b)$p.value

  # if checking for normality do a wilcox test instead of t.test
  res = 'if'(chk_norm & (norm_pA < norm_a | norm_pB < norm_a),
             wilcox.test(a, b)$p.value,
             t.test(a, b, var.equal = T)$p.value)
  return(res)
}


# Run sim for the 4 scenarios
res = list()
res[[1]] = replicate(1e4, sim(chk_norm = T, rem_out  = T))
res[[2]] = replicate(1e4, sim(chk_norm = F, rem_out  = T))
res[[3]] = replicate(1e4, sim(chk_norm = T, rem_out  = F))
res[[4]] = replicate(1e4, sim(chk_norm = F, rem_out  = F))

# Inspect % p under 0.05 for the 4 scenarios
sapply(res, function(x) mean(x < 0.05))

# Look at null distributions under the 4 scenarios
par(mfrow = c(2, 2))
for(i in 1:length(res)){
  hist(res[[i]], col = "Light Blue",xlab = "p", main = paste("Scenario", i))
}