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$?


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) {
  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)
  • $\begingroup$ To get a specified prevalence, you're going to need to modulate both the distributions of your X variables (here binary), and your betas (here group differences). With 1k X's, it's going to be messy. You also need to consider how you want your X's to be related to each other. $\endgroup$ Dec 3, 2019 at 20:25
  • $\begingroup$ @gung-ReinstateMonica - I meant to add (now added) that the binary variables are all independent. $\endgroup$ Dec 3, 2019 at 20:27
  • 1
    $\begingroup$ Fair enough. You basically have 3 moving parts, but 2 degrees of freedom. Start by deciding the rate you want for the intercept, then (say) the prevalences of each X. Then the required betas will be fixed. Or specify the betas & the prevalences will be fixed. Bear in mind that with many X's, there will be many combinations of beta values that will sum to the required total. $\endgroup$ Dec 3, 2019 at 20:31
  • $\begingroup$ @gung-ReinstateMonica - Thanks! That helped me to clear up my thinking and get a working simulation. $\endgroup$ Dec 3, 2019 at 22:32
  • $\begingroup$ You're welcome, @stats134711. If you want, I can try to write up something perfunctory tomorrow that summarizes my comments. $\endgroup$ Dec 4, 2019 at 3:44

1 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).


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.