Simulating population-level logistic regression model with pre-specified prevalence

Question

I'm interested in simulating the following prospective, population-based model for binary outcome $Y_i$, and independent subjects $i=1,\dots,N$: $$ \Pr(Y_i=1\mid X_i,G_i)=\frac{1}{1+\exp(-(\alpha_0 + 0.5 X_i + G_i\beta))} $$ where $X_i\sim N(0,1)$ is a continuous covariate with effect 0.5 (log OR scale), $G_i$ is an $m\times 1$ vector of binary covariates, with corresponding effects $\beta$. The value of $m$ can be large and is usuall around 1000. Each element of $G_i$ is generated independently with some probability $\pi_j$, $j=1,\dots,m$. That is $\mathbb{E}G_{ij}=\pi_j$. $\alpha_0$ is an intercept related to the prevalence, $p_0=\Pr(Y_i=1)$ of outcome in the population.

Question When I run the following model below, I set $\alpha_0=\log p_0/(1-p_0)-\sum_{j=1}^m\pi_{j}\beta_j$ in the probability model above. The idea of subtracting out the second term is so that on average it can cancel out with $G_i\beta$, leaving just the log odds (first term).

However, I find that my estimated prevalence $\hat{\Pr}(Y_i=1)=\sum_{i=1}^N I(Y=1)/N$=0.27948 is not even close to the value of $p_0=0.01$. This happens because $\sum_{j=1}^m\pi_{j}\beta_j$ is a non-negligible term. How can I set my intercept to get the right estimate?

The probability model I have is a conditional probability and $p_0=\Pr(Y_i=1)$. Taking the expectation with respect to $X$ and $G$ would yield the marginal, but it can get quite messy. Is there a slightly adhoc way just for simulation purposes to get close to the desired $p_0$?

EDIT

Using gung's comments, I've modified the code below so that I fix $\pi_j$'s and $\beta$'s first. Then $\alpha_0$ is searched until the estimate $\hat{\Pr}(Y_i=1)\approx p_0$.

#Population size
N = 100000

#Prevalence of disease
p0 = 0.01

#Continuous covariate
X <- rnorm(N,0,1)

#Number of binary covariates
m = 1000

#Set pi_j to be the same for each binary covariate
pi_j = rep(0.01,m)

#Generate binary covariates for N subjects
#See below for genSpMat definition
G <- genSpMat(nrows=N,ncols=m,col_probs=pi_j)

#Generate betas for each binary covariate
betas <- runif(m,1,4)

#Generate probability model
#Guess alpha0 so that E(pr)~0.01
#Max alpha0 is log(p0/(1-p0))

alpha0 <- log(p0/(1-p0))

pr <- rep(1,N)
#Adjust intercept until population prevalence is reached (to a certain tolerance)
while(abs(mean(pr)-p0) > 2*10^(round(log10(p0))-1)){
  alpha0 <- alpha0 - 1
  pr <- plogis(alpha0 + 0.5*X + as.vector(G%*%betas))
}

Y <- rbinom(N,1,pr) #Outcome


#Generate sparse matrix of 0/1s
genSpMat <- function(nrows, ncols, col_probs) {
  require(Matrix)
  r <- lapply(1:ncols, function(x) {
    p <- col_probs[x]
    i <- sample.int(2L, size = nrows, replace = T, prob = c(1 - p, p))
    which(i == 2L)
  })
  rl <- lengths(r)
  nc <- rep(1:ncols, times = rl) # col indexes
  nr <- unlist(r) # row index
  ddims <- c(nrows, ncols)
  #ngTMatrix output
  sparseMatrix(i = nr, j = nc, dims = ddims,giveCsparse = FALSE)
}

Answer 1

If you want the resulting prevalence to have a specific (population) value, you will need to modulate both the distributions of your $X$ variables, and your betas. You will also need to decide what rate you want for the intercept, and you need to take into account the pattern of correlations amongst your $X$ variables. In essence, you have several moving parts, but fewer degrees of freedom (for a simpler, analogous example, it may help to read my answer here: How to simulate censored data). Thus, you will need to stipulate some of these, and solve for the rest. Note also that given some preferred values for the first few of these, it may not be possible to find values for the last that will get you to the result you want (for examples of this, see the prior answer or here: How to simulate data that satisfy specific constraints such as having specific mean and standard deviation?). Lastly, note that with multiple $X$'s, there will be an infinite number of possible solutions, for example, possible betas for $X_1$ and $X_2$ could be:
\begin{align} \beta_1\qquad &\beta_2 \\ \hline \\[-10pt] .5\qquad &.2 \\ .4\qquad &.3 \\ .3\qquad &.4 \\ .2\qquad &.5 \\ \!\vdots\qquad &\;\,\vdots \end{align}

Your situation seems a little bit simpler, in that your $X$'s are all binary, and they are all independent. You will still need to stipulate values for some of these moving parts, and solve for the last ones. If I were doing this, I would start by specifying the rate you want for the intercept, and the prevalences of each $X$, then I would generate betas subject to some constraints (e.g., all positive, and must sum to a given limit).