# using factors (resulting in mean(Group)=0.5, mean(Sex)=0.5)
sim = function(n=40, mu=0.6343) {
  a  = rnorm(n,   0.0, 1)
  b1 = rnorm(n/2, mu - mu/4, 1)
  b2 = rnorm(n/2, mu + mu/4, 1) # interaction is mu/2

  dat = data.frame(Group = c(rep("C", n), rep("T", n)),
                   Sex   = rep(c(rep("F", n/2), rep("M", n/2)), 2),
                   Val   = c(a, b1, b2))
  fit   = lm(Val ~ Group*Sex, dat)
  return(summary(fit)$coefficients)
}

# mean(Sex)=0, mean(Group)=0
sim2 = function(n=40, mu=0.6343) {
  a  = rnorm(n,   0.0, 1)
  b1 = rnorm(n/2, mu - mu/4, 1)
  b2 = rnorm(n/2, mu + mu/4, 1) # interaction is mu/2

  dat = data.frame(Group = c(rep(-0.5, n), rep(0.5, n)),
                   Sex   = rep(c(rep(-0.5, n/2), rep(0.5, n/2)), 2),
                   Val   = c(a, b1, b2))
  fit   = lm(Val ~ Group*Sex, dat)
  return(summary(fit)$coefficients)
}

n = 40
d = power.t.test(n=n, power=0.8)$delta

res=sapply(1:1000, function(x) sim(n, d)) # variant 1
#res=sapply(1:1000, function(x) sim2(n, d)) # variant 2

mean(res[14,]<0.05) # real power main effect
mean(res[16,]<0.05) # real power interaction
data.frame(sd.of.coef.sample=sapply(1:4, function(x) sd(res[x,])),
           mean.coef.se=sapply(5:8, function(x) mean(res[x,])),
           q25.coef.se=sapply(5:8, function(x) quantile(res[x,], 0.25)),
           q75.coef.se=sapply(5:8, function(x) quantile(res[x,], 0.75)),
           coeff=sapply(1:4, function(x) mean(res[x,])),
           z.coeff=sapply(1:4, function(x) mean(res[x,]/mean(res[x+4,]))))

mean(res[2, res[14,]<0.05]) # mean of significant main coefs
mean(res[4, res[16,]<0.05]) # mean of significant interaction coefs

mean(res[2, res[14,]<0.05]) / d     # inflation main coef
mean(res[4, res[16,]<0.05]) / (d/2) # inflation interaction coef