# Generate correlated binary variable from correlated normal data

I have a correlated normal data

y <- -1 + 3*x1- 2*x2 + e


where, e is generated from multivariate normal distribution and thus correlated.

To get a correlated binary data I followed these steps

1. scale to get a marginally standard Normal distribution.
2. transform by their cumulative distribution function to get a uniform distribution
3. use the inverse transform method to get binomial responses.

st.y <- (y-mean(y))/sd(y)

erfc <- function(x) 2 * pnorm(x * sqrt(2), lower = FALSE) #complementary error fuction
u <- .5*erfc(-st.y/sqrt(2))

z <- ifelse(u<=.5, 0,1)


But if I fit gee from this data the pre-specified parameter of intercept and slope is changed. Because in glm, logit(z) = $\beta_0$ + $\beta_1$x1 + $\beta_2$x2.

How can I get a correlated binary variable satisfying logit(z) = -1+3x1-2x2?