How to set up a DGP for Monte Carlo simulation with non-independent regressions (correlated errors) I want to set up a data generating process for two different estimations. The idea is to show how bias is introduced when the models are not properly specified. The first model should be a logit/probit estimation and the second one a form of a multinomial regression. At first I thought that this is straight forward but I seem to be unable to correlate the two regressions according to my liking.
If, say, the regressions were both linear, then I would proceed like this:
set.seed(12)
x1 <- runif(100) * 4 - 2

beta0 <- c(0, 2)
beta1 <- c(1, -1)

#setting up the correlated errors:
vmat <- matrix(c(1,0.7, 0.7, 1), nrow=2)
v <- mvtnorm::rmvnorm(100, c(0,0), vmat)

e1 <- v[,1]
e2 <- v[,2]

lin.pred1 <- beta0[1] + beta1[1]*x1 + e1
lin.pred2 <- beta0[2] + beta1[2]*x1 + e2

fit1 <- lm(lin.pred1~x1)  
fit2 <- lm(lin.pred2~x1)  

cor(residuals(fit1), residuals(fit2))  
[1] 0.6729155

However, how does one correlate errors for a logit/probit with another regression?
Best,
Tom
 A: I can think of two ways to do this, and there may be more.
The first way involves including an unobserved predictor in both models that is not correlated with other variables (so it does not induce confounding). This is a very simple method that accomplishes the task of inducing a correlation between the two outcomes, but it is not straightforward to control what that correlation is (which depends on how you measure it anyway) or to compute the coefficients of a model that marginalizes over the additional predictor, as the true coefficients of the other predictors will be different from what you set them to be. You would have to use simulation to figure out the true values of the coefficients in a model that ignores the additional predictor.
The second way involves setting a probit regression model using the latent variable formulation. You can generate the latent errors using a multivariate normal distribution just as you did in the linear case and then supplying a threshold to make the outcomes binary. This is the same data-generating process used for generating data for categorical confirmatory factor analysis, so there may be literature there that could be helpful.
