# Simulate study of treatment vs control w/ half male subjects
sim = function(n, show_dat = FALSE, show_fit = FALSE, sex_main = F, intercept = F){
  a  = rnorm(n,   0.0, 1)
  b1 = rnorm(n/2, 0.75, 1)
  b2 = rnorm(n/2, 1.25, 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))

  if(sex_main){
    fit = 'if'(intercept, lm(Val ~ 0 + Group*Sex, dat),         lm(Val ~ Group*Sex, dat))
  }else{
    fit = 'if'(intercept, lm(Val ~ 0 + Group + Group:Sex, dat), lm(Val ~ Group + 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 Simulations
par(mfrow = c(2, 2))
nvals = seq(4, 100, by = 2)
nsim  = 1e3
out   = list()
mat   = expand.grid(c(TRUE, FALSE), c(TRUE, FALSE))
for(i in 1:nrow(mat)){
  res = out[[i]] = sapply(nvals, function(n)
                          rowMeans(replicate(nsim, sim(n, sex_main = mat[i, 1], intercept = mat[i, 2])) < 0.05))

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

  # Simulation Results
  main = paste0(c("sex_main:", "intercept:"), mat[i,])
  plot(nvals, res[1, ], type = "n", ylim = range(res), xlab = "n", ylab = "% Sig", main = main, 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/2, .5), lwd = 2, lty = 2)

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


  res[, nvals == n80]
  exactPwr(n80/2, .5)
}
cbind(exactPwr(nvals, 1), exactPwr(nvals/2, .5))