# Simulate data for power analysis of logistic regression model - include covariance variance of variables?

I've tried to simulate data for a power analysis of a logistic regression. The results of the power analysis look reasonable: power=90% for a sample of 6000 persons. But I feel that the analysis lacks something. So, my question is: when generating the data should I include something about how the variables are correlated, or their covariance, other than just defining their linear relationship as I have done in the example below, and if so where do I write that into the code?

I know other questions look like this but I'm not confident that they answer this question.

library(plyr) # functions
## Define Function
simfunktion <- function() {
# Number in each sample
antal <- 6000
beta0 <- log(0.16) # logit in reference group
beta1 <- log(1.1)  # logit given smoking
beta2 <- log(1.1)  # logit given SNP(genevariation)
beta3 <- log(1.2)  # logit for interactioncoefficient for SNP*rygning
## Smoking variable, with probabilities defined according to empirical studies.
smoking  <- sample(x = 0:2, size = antal, replace = TRUE, prob = c(40,25,40))
## SNP variables with probabilities defined according to empirical studies
SNP      <- sample(x = 0:2, size = antal, replace = TRUE, prob = c(40,40,20))
## calculated probabilites given the model:
pi.x     <- exp( beta0 + beta1*smoking + beta2*SNP + beta3*smoking*SNP) /
( 1 + exp(beta0 + beta1*smoking + beta2*SNP + beta3*smoking*SNP) )
## binoial events given the probabilities:
sim.y    <- rbinom( n = antal, size = 1, prob = pi.x)
sim.data <- data.frame(sim.y, smoking, SNP)
#################### p-value of the interaction is extracted:
## the model is run:
glm1     <- glm( data = sim.data, formula = sim.y ~ smoking + SNP + smoking:SNP,
family=binomial )
## p-value of the interactionterm is extracted:
summary(    glm( data = sim.data, formula = sim.y ~ smoking + SNP + smoking:SNP,
family=binomial ))\$coef[4,4]
}
pvalue     <- as.vector(replicate( 100 , simfunktion()))
mean(pvalue < 0.05)

Let me throw out some thoughts, and we'll see if something helps you.

Some preliminaries:

• I find it a bit odd that you are defining your true betas as the log of some number; is that because you are using reported odds ratios? (If so, this is perfectly appropriate.)
• It's important to realize when doing power analyses based on effects reported in the literature that the results are optimistic. I discuss that here: Desired effect size vs. expected effect size; you may also want to read this thread: Logistic regression model manipulation.
• I notice that you are simulating the expected distribution of your covariates. That's not typically done; in general, we assume that our covariates are a set of known constants. However, if you will be doing observational (e.g., epidemiological) research, these can well vary and this strategy is appropriate.
• Your covariates have the values 0, 1, 2. Are these levels of a factor, or are they equal interval? I ask because they look like the former, but are treated as the latter in the data generating process.