# Simulate study of treatment vs control w/ half male subjects
sim = function(n, show_dat = FALSE, show_fit = FALSE){
  a  = rnorm(n,   0.0, 1)
  b1 = rnorm(n/2, 0.5, 1)
  b2 = rnorm(n/2, 1.5, 1)

  dat = data.frame(Group = c(rep("C", n), rep("T", n)),
                   Sex   = rep(c(rep("M", n/2), rep("F", n/2)), 2),
                   Val   = c(a, b1, b2))

  fit   = lm(Val ~ Group*Sex, dat)
  pvals = summary(fit)$coefficients[, 4]

  if(show_dat) print(dat)
  if(show_fit) print(summary(fit))
  return(pvals)
}

# Analytical power calculation by effect and sample size
exactPwr = function(nvals, delta){
  sapply(nvals, function(n)
         power.t.test(n = n, delta = delta, sd = 1, sig.level = .05, strict = T)$power)
}


# Run Simulation
nvals = seq(4, 100, by = 2)
res   = sapply(nvals, function(n) rowMeans(replicate(1e4, sim(n)) < 0.05))

# Get Sample size nearest to 80% power
n80  = nvals[which.min(abs(res["GroupT",] - 0.8))]


# Simulation Results
plot(nvals, res[1, ], type = "n", ylim = range(res), xlab = "n", ylab = "% Sig", panel.first = grid())
for(i in 1:nrow(res)){
  lines(nvals, res[i,], col = rainbow(4)[i],  lwd = 3)
}

# Add analytical results
lines(nvals, exactPwr(nvals, 1),  lwd = 2, lty = 2)
lines(nvals, exactPwr(nvals, .5), lwd = 2, lty = 2)

abline(h = 0.8, v = n80, lty = 2)
legend("right", c(rownames(res), "Analytic Match"),
       col = c(rainbow(4), "Black"), lwd = 3, lty = c(rep(1, 4), 2), bty = "n")


res[, nvals == n80]
exactPwr(nvals, .5)[nvals == n80]