# Simulate event rate in logistic regression - Finding the intercept

I want to simulate a logistic regression (using a set of continous and binary confounders with known odds ratios) with a specified probability outcome (e.g. event rate = 0.2).

This is actually a somewhat follow-up question to: Simulate rare event data using logistic regression with correlated covariates in R. (Intuitively, I would have posted a comment on this but unfortunately, I have not enough reputation to do so. Following a discussion in meta on this issue, there is an advise to ask a new question.)

So, to my understanding, simulating the probability rate when the regression coefficients are known must be solved by finding the correct intercept. When I was trying on this before reading that thread, I had considered an alternate, rather simple approach:

$log(P/(1-P ) = \beta_0 + \beta_1*X_1 + ... + \beta_n*X_n$

and therefore,

$\beta_0 = log(P/1-P) - \beta_1*X_1 - ... - \beta_n*X_n$

Q1) Obviously, this approach has a major flaw (see R Code attached). Can anyone explain me my errors in reasoning?

whuber has proposed an interesting approach (unfortunately, I have some problems of comprehense) to set σ=1 and solve numerically for $β_0$:

$\mathbb{E}(\Pr(Y=1|X)) = \frac{1}{\sqrt{2\pi}}\int_\mathbb{R} \left(\frac{\exp(-z^2/2)}{1 + \exp(-\beta_0 - \sigma z)}\right)dz.$

He also presented a plot of Pr(Y=1) as a function of $β_0$ with σ=1.

I would be most grateful if someone could explain to me:

Q2.) why you obtain the expectation by integrating against the density $\exp(-z^2/2)/\sqrt{2\pi}$

Q3.) how to solve this equation for $\beta_0$ and

Q4.) how this plot was generated.

# set.seed(1)

# Sample size
n <- 10^4

# Confounders (2 contin. + 2 binary vars)
cont_var <- replicate(2, rnorm(n))
bin_var <- replicate(2, rbinom(n, size = 1, p = 0.5))

df <- data.frame(cont_var, bin_var)

# Set Odds ratios (0.1 1.4 2.7 4.0)
beta_var <- seq(0.1, 4.0, length.out = 4)

# Multiply betas with confounders
i <- 1 : ncol(df)
x <- mapply(*, log(beta_var[i]), df[i])

# Let's assume the Event rate = 0.2
event_rate <- 0.2

# Find beta_0
beta_0 <- log(event_rate/(1 - event_rate)) - rowSums(x)
beta_0 <- mean(beta_0)

# Generate Outcome
z <- cbind(beta_0, x)
z <- rowSums(z)
prob <- 1/(1 + exp(-z))    # Inverse logit
outcome <- rbinom(n, size=1, prob=prob)

# Check OR
d1 <- data.frame(df, outcome)
model <- glm(outcome ~ ., data = d1, family = "binomial")
exp(coef(model)) # correct!

# Check if Event_rate is correct
sum(outcome) / length(outcome) # not correct!