Following on from the post by Stephan Kolassa above (I can't add this as a comment), I have some alternative code for a simulation.  This uses the same basic structure, but is exploded a bit more.  It also is based on the code by [Kleinman and Horton][1] to simulate the logistic regression.  

nn is the number in the sample.  The covariate should be continuously normally distributed, and standardized to mean 0 and sd 1.  We use rnorm(nn) to generate this.  We select an odds ratio and store it in odds.ratio.  We also pick a number for the intercept.  Choice of this number governs what proportion of the sample experience the "event" (e.g. 0.1, 0.4, 0.5).  You have to play around with this number until you get the right proportion.  The following code gives you a proportion of 0.1 with a sample size of 950 and an OR of 1.5:

    nn <- 950
    runs <- 10000
    intercept <- log(9)
    odds.ratio <- 1.5
    beta <- log(odds.ratio)
    proportion  <-  replicate(
                  n = runs,
                  expr = {
                      xtest <- rnorm(nn)
                      linpred <- intercept + (xtest * beta)
                      prob <- exp(linpred)/(1 + exp(linpred))
                      runis <- runif(length(xtest),0,1)
                      ytest <- ifelse(runis < prob,1,0)
                      prop <- length(which(ytest <= 0.5))/length(ytest)
                      }
                )
    summary(proportion)

 
summary(proportion) confirms that the proportion is ~ 0.1

Then using the same variables, the power is calculated over 10000 runs:

    result <-  replicate(
                  n = runs,
                  expr = {
                      xtest <- rnorm(nn)
                      linpred <- intercept + (xtest * beta)
                      prob <- exp(linpred)/(1 + exp(linpred))
                      runis <- runif(length(xtest),0,1)
                      ytest <- ifelse(runis < prob,1,0)
                      summary(model <- glm(ytest ~ xtest,  family = "binomial"))$coefficients[2,4] < .05
                      }
                )
    print(sum(result)/runs)


I think that this code is correct - I checked it against the examples given in Hsieh, 1998 (table 2), and it seems to agree with the three examples given there.  I'd love any feedback on the code, and examples of where it breaks down.  


  [1]: http://sas-and-r.blogspot.com/2009/06/example-72-simulate-data-from-logistic.html