Logistic regression simulation with respect to event occurrence (prevalence) I am trying to simulate logistic regression data, but under the constraints of prevalence.
$$\text{logit}(y_i) = \beta_0 + \beta_1 X_1 + \beta_2X_2$$
For example, I want to create a dataset that has 10% of $y = 1$, and 90% $y = 0$.
I can think of two approaches:

*

*search a grid of $\beta_0$, that produce the closest prevalence to 10%

*post hoc resampling. Stage 1: fix all $\beta$s to produce $y$. Stage2: sampling from y so that prevalence to be 10%.

Here is the code for approach 2:
# fix betas
betas = c(-2, 1, -0.5)

# X could be anything, I just use binomial here
# dimension = 10,000 * 3
X = cbind(1, matrix(rbinom(20000, size = 2, prob = 0.2), ncol = 2))

# generate probability according to sigmoid function
probs = 1/(1 + exp(- X %*% betas))

# generate y
y = apply(probs, 1, function(prob){rbinom(1, size = 1, prob = prob)})
  
# resample y, so that prevalence can be constrained
# dimension = 1,000 * 3
case = sample(which(y == 1), size = 100)
control = sample(which(y == 0), size = 900)
  
resampleX = X[c(case, control), ]
resampley = y[c(case, control)]
# run glm 
glm(resampley ~ 0 + resampleX, family = binomial)$coefficients

 
 
The question is: which approach is valid?
 A: You have an array of explanatory variables $(x_1, x_2, \ldots, x_n)$ ($n=20000$) and a model that assigns a probability to each $x_i.$  You seek a subarray of these variables that has a mean probability of $p=0.10.$  Clearly you don't need that mean to be exactly $0.10$--that might be impossible to achieve and would be far too precise a demand in any actual situation.  So, how can you come reasonably close to $0.10$?
The most general solution selects $x_i$ independently with a chance $q_i.$  (Set $q_i=1$ to specify a definite selection and $q_i=0$ to omit $x_i.$)  The expected mean probability then is
$$\mu(\mathbf q) = \sum_{i=1}^n p_i q_i$$
and the expected number of $x_i$ chosen is $\sum q_i.$
You could, for instance, select probabilities from smallest to largest until reaching the desired proportion.  This would maximize the number of retained observations. Such methods, though, leave us with the uneasy feeling that your selected data might not represent the phenomenon you are trying to study: you are choosing only the kinds of observations that lead to probabilities near your target and are excluding observations with the highest probabilities.
One relatively simple way to alleviate this concern is to find a number $k$ for which $q_i \propto \exp(-k p_i)$ works and, in order to obtain as many $x_i$ as possible, to maximize the constant of proportionality to make the largest $q_i$ equal to $1.$  This probability sampling solution enables you to relate properties of your sample to properties of the original $X$, if you wish; and it offers a chance to select observations yielding the whole spectrum of response values.
Here's an R implementation. It begins after you have specified the array X and computed the probabilities probs according to your model.  It finds $k$, uses that to compute the $q_i,$ and then stores a random selection in the array of indexes i.  The dataset of selected explanatory variables can then be obtained as X[i, ]
f <- Vectorize(function(k, p, target=0.10) sum(p * exp(-k*p)) / sum(exp(-k*p)) - target, "k")
uniroot(f, c(0, log(10^308)), p = probs)$root -> k
q <- exp(-k * probs) 
i <- runif(length(probs)) * max(q) < q

For instance, after running set.seed(17), the mean value of the selected probabilities was $0.1003$ and $n^\prime=2973$ rows of X were selected.  This is not too bad, given that over three quarters of your original $p_i$ exceed the target of $0.1,$ and the most you can select is $5583$ in this example.
This approach preserves much of the variation in the probabilities, too, as you can see by comparing its histogram (middle) to the solution with the $5583$ smallest probabilities (right):

Of course, the dataset you ultimately generate will depend randomly on these adjusted probabilities.  Without loss of generality, let the selected observations have indexes $1,2,\ldots, n^\prime.$  The expectation of the mean response will be equal to $\mu(\mathbf q)$ (very nearly your target value) while its variance will be
$$\operatorname{Var}\left(\frac{1}{n^\prime}\sum_{i=1}^{n^\prime} Y_i\right) =\frac{1}{{n^\prime}^2} \sum_{i=1}^{n^\prime}p_i(1-p_i) \le \frac{\mu(\mathbf q)}{{n^\prime}}.$$
Thus, when this procedure selects an adequately large subset of $n^\prime$ values, the proportion of positive responses in the dataset itself will therefore be very close to your target value, differing from it due to its random construction:
y <- runif(sum(i)) <= probs[i]
mean(y)


0.09496

A: Its going to be hard to simulate a with an exact proportion of 1s, but if you can get pretty close if you simulate a lot of data.
Key thing to realize is that the intercept is the log odds of the outcome when the covariates are mean centred or have 0 as the reference.  So let
$$ \beta_0 = \log(0.1/0.9) $$
and simulate a bunch of data.  Shown below is some code to do this which gets fairly close to your requirement.
set.seed(0)
N <- 100000
x1 <- rbinom(N, 1, 0.5)
x2 <- rbinom(N, 1, 0.5)
intercept <- qlogis(0.1) 
b1 <- 1
b2 <- -2
beta <- c(intercept, b1, b2)
X <- model.matrix(~x1 + x2)
eta <- X %*% beta
y <- rbinom(N,1, plogis(eta))

mean(y)
[1] 0.09622

```

