I looked into resources that graciously provided code for

-simulating logistic regression, (e.g. http://en.wikibooks.org/wiki/R_Programming/Binomial_Models)

-simulating correlated multivariate normals (e.g. R package "MSBVAR")

I would like to simulate a setting of n.obs = 100,000 with 5% probability of the binary outcome Y being 1 and the correlation of 4 normally distributed continuous covariates V1 V2 V3 and V4 being: Pearson's product moment correlation coefficients

My question is: how could I determine the value of the coefficients for the covariates?

Any help is welcome. Thank you for your input!

Edit 03/17/2015 for Question on the value of beta's after the βΣβ′ step: Is there a way to obtain larger beta values after the βΣβ′ step while having them equal to each other?enter image description here


1 Answer 1


Let the variables be $X_1, \ldots, X_4$. Due to the flexibility in choosing coefficients in the regression model, we will lose no generality by assuming they all have zero means. Write the covariance matrix as $\Sigma$. Let the logistic regression coefficients be $\beta_0, \beta_1, \ldots, \beta_4$. This means the model is

$$\text{logit}(\Pr(Y=1)) = \beta_0 + \beta_1 X_1 + \cdots + \beta_4 X_4.$$

The right hand side, as an affine combination $X$ of components of a multivariate Normal vector, thereby has a Normal distribution. All we need to know are its mean, which equals $\beta_0$, and its variance, which is

$$\text{Var}(\beta_0 + \beta_1 X_1 + \cdots + \beta_4 X_4) = \beta\Sigma \beta^\prime,$$

with $\beta=(\beta_1, \beta_2, \beta_3, \beta_4)$. For future reference let's call the right hand side $\sigma^2$.

The unconditional probability $\Pr(Y=1)$ is the expectation of $\Pr(Y=1|X)$. To find this, re-express the preceding result in terms of the probability directly:

$$\Pr(Y=1) = \frac{1}{1 + \exp(-\beta_0 - \sigma Z)}$$

where $Z$ has a standard Normal distribution. Obtain the expectation by integrating against the density $\exp(-z^2/2)/\sqrt{2\pi}$:

$$\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.$$

This has to be solved numerically for $\sigma$ and $\beta_0$, giving a one-dimensional manifold of solutions. A simple way to obtain a solution is to set $\sigma=1$ and solve numerically for $\beta_0$. Alternatively, choose $\beta_1, \ldots, \beta_4$ any way you wish, thereby determining $\sigma$, and then solve for $\beta_0$. (The advantage of setting $\sigma$ to a standard value, such as unity, is that you can tabulate values of $\beta_0$ for a wide range of possible values of $\Pr(Y=1)$, once and for all, obviating any need to incorporate the numerical integration and root-finding code in every application.)

Figure--plot of the expectation vs. beta

This figure plots $\Pr(Y=1)$ as a function of $\beta_0$, with $\sigma=1$. Because it rises continuously and monotonically from $0$ (when $\beta_0\to-\infty$) to $1$ (when $\beta_0\to+\infty$), every possible value of $\Pr(Y=1)$ can be realized by a unique corresponding value of $\beta_0$.

The solution for $\Pr(Y=1)=0.05$ is approximately $\beta_0=-3.37$, where the graph (in the figure) attains a height of $0.05$ (as shown with a dotted line).

This provides infinitely many solutions, because $\beta_1, \ldots, \beta_4$ can be anything you like: selecting them determines $\sigma$, from which $\beta_0$ can be computed.

The correctness and accuracy of this approach are supported by an R-based simulation. It generates $\beta_1,\ldots,\beta_4$ randomly and uses the previously computed value of $\beta_0$ to construct a dataset of 100,000 records. At the end it displays the correlation matrix of the independent variables (to check it has the desired coefficients) and it outputs the mean of the dependent variable to see whether it exhibits $5\%$ positive results overall (at least up to chance variation).

Repeated runs of this simulation (which thereby involve all new variable values as well as a new model in each iteration) consistently produce proportions of positive results between $4.8\%$ and $5.2\%$. In 100 iterations (starting with a seed of $17$), the average proportion was $4.99\%$.

library(MASS) # exports mvrnorm()
# Describe the independent variables.
mu <- c(0,0,0,0)
Sigma <- cbind(c(10,6,5,6), c(6,10,5,3.5), c(5,5,10,3), c(6,3.5,3,10))/10
# Simulate a dataset in which Pr(Y=1) is 5%.
beta.0 <- -3.37154                  # Corresponds to 5%.
n.obs <- 1e5
x <- mvrnorm(n.obs, mu, Sigma)
beta <- rnorm(4)                    # Can be anything (nonzero)!
sigma2 <- beta %*% Sigma %*% beta  
beta <- beta / sqrt(sigma2)         # Assures sum of squares is unity
y <- runif(n.obs) < 1 / (1 + exp(-beta.0 - x %*% beta))
# Confirm that the independent variables have the desired correlation
# and the dependent variable has the desired proportion of true responses.
round(cor(x), 2)
  • $\begingroup$ Thank you whuber! I just tried the code a few times and it runs beautifully. Is there a way to constrain the final beta to be within a certain range? (like all beta's being positive?) $\endgroup$
    – HueSX
    Mar 10, 2015 at 15:55
  • 1
    $\begingroup$ The only constraint is that $\beta_0\lt 0$ for $\Pr(Y=1) \lt 1/2$ and $\beta_0\gt 0$ for $\Pr(Y=1) \gt 1/2$. This holds regardless of the value of $\sigma$. $\endgroup$
    – whuber
    Mar 10, 2015 at 17:30
  • $\begingroup$ Hi Mr whuber May I ask a follow-up question about your code? I tried different combinations of initial beta values, and it seems that the the solution gets as low as around 0.321 regardless of initial values. Could you please give me some pointers on why this is the case and if I wish to yield beta's of higher values (say around 1) what to do? Thank you in advance! $\endgroup$
    – HueSX
    Mar 17, 2015 at 23:20
  • $\begingroup$ I'm afraid I don't understand your question. What exactly do you mean by "beta values" (beta or beta.0?) and what do you mean by "the solution" (the simulation, the correlation coefficient, the mean of y, something else)? $\endgroup$
    – whuber
    Mar 17, 2015 at 23:23
  • $\begingroup$ I tried for example starting with beta's all equal to 1 or all equal to 4 or other combinations (did not change beta.0 in your code), then after the βΣβ′ step, the β's equal 0.321-ish. Question please see above (incorporated into original post) $\endgroup$
    – HueSX
    Mar 18, 2015 at 0:10

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.