Simulation of logistic regression power analysis - designed experiments This question is in response to an answer given by @Greg Snow in regards to a question I asked concerning power analysis with logistic regression and SAS Proc GLMPOWER.
If I am designing an experiment and will analze the results in a factorial logistic regression, how can I use simulation ( and here ) to conduct a power analysis?
Here is a simple example where there are two variables, the first takes on three possible values {0.03, 0.06, 0.09} and the second is a dummy indicator {0,1}. For each we estimate the response rate for each combination (# of responders / number of people marketed to). Further, we wish to have 3 times as many of the first combination of factors as the others (which can be considered equal) because this first combination is our tried and true version. This is a setup like given in the SAS course mentioned in the linked question.

The model that will be used to analyze the results will be a logistic regression, with main effects and interaction (response is 0 or 1). 
mod <- glm(response ~ Var1 + Var2 + I(Var1*Var2))

How can I simulate a data set to use with this model to conduct a power analysis?
When I run this through SAS Proc GLMPOWER (using STDDEV =0.05486016
 which corresponds to sqrt(p(1-p)) where p is the weighted average of the shown response rates):  
data exemplar;
  input Var1 $ Var2 $ response weight;
  datalines;
    3 0 0.0025  3
    3 1 0.00395 1
    6 0 0.003   1
    6 1 0.0042  1
    9 0 0.0035  1
    9 1 0.002   1;
run;

proc glmpower data=exemplar;
  weight weight;
  class Var1 Var2;
  model response = Var1 | Var2;
  power
    power=0.8
    ntotal=.
    stddev=0.05486016;
run;

Note: GLMPOWER only will use class (nominal) variables so 3, 6, 9 above are treated as characters and could have been low, mid and high or any other three strings. When the real analysis is conducted, Var1 will be used a numeric (and we will include a polynomial term Var1*Var1) to account for any curvature.
The output from SAS is 

So we see that we need 762,112 as our sample size (Var2 main effect is the hardest to estimate) with power equal to 0.80 and alpha equal to 0.05. We would allocate these so that 3 times as many were the baseline combination (i.e. 0.375 * 762112) and the remainder just fall equally into the other 5 combinations.
 A: @Gung's answer is great for understanding.  Here is the approach that I would use:
mydat <- data.frame( v1 = rep( c(3,6,9), each=2 ),
    v2 = rep( 0:1, 3 ), 
    resp=c(0.0025, 0.00395, 0.003, 0.0042, 0.0035, 0.002) )

fit0 <- glm( resp ~ poly(v1, 2, raw=TRUE)*v2, data=mydat,
    weight=rep(100000,6), family=binomial)
b0 <- coef(fit0)


simfunc <- function( beta=b0, n=10000 ) {
    w <- sample(1:6, n, replace=TRUE, prob=c(3, rep(1,5)))
    mydat2 <- mydat[w, 1:2]
    eta <- with(mydat2,  cbind( 1, v1, 
                v1^2, v2,
                v1*v2,
                v1^2*v2 ) %*% beta )
    p <- exp(eta)/(1+exp(eta))
    mydat2$resp <- rbinom(n, 1, p)

    fit1 <- glm( resp ~ poly(v1, 2)*v2, data=mydat2,
        family=binomial)
    fit2 <- update(fit1, .~ poly(v1,2) )
    anova(fit1,fit2, test='Chisq')[2,5]
}

out <- replicate(100, simfunc(b0, 10000))
mean( out <= 0.05 )
hist(out)
abline(v=0.05, col='lightgrey')

This function tests the overall effect of v2, the models can be changed to look at other types of tests.  I like writing it as a function so that when I want to test something different I can just change the function arguments.  You could also add a progress bar, or use the parallel package to speed things up.
Here I just did 100 replications, I usually start around that level to find the approximate sample size, then up the itterations when I am in the right ball park (no need to waste the time on 10,000 iterations when you have 20% power).
