How to account for Covariances in a Monte Carlo simulation?

Question

I am conducting a regression analysis of Z on X and Y. I am interesting in interpreting the marginal effect. The model includes quadratic terms for both X and Y as well as an interaction term X*Y. If I look at just the marginal effect of X on Z (the partial derivative of Z wrt X) I get the equation:

ME <- function(x, y){0.187 - 0.000472*x - 0.000252*y}

I am interested in the point at which this ME becomes zero and hence the partial effect becomes negative. I am calling this a "tipping point". So to calculate the tipping point, let the ME equation follow the form: 0 = a + bx +cy. Solving this for x gives the tipping point im interested in. Rather than taking this at face value, I would like to construct a distribution using a Monte Carlo simulation. I created the following r code:

set.seed(123)
a = rnorm(1000, mean=0.187, sd=0.016711349)
b = rnorm(1000, mean=-0.000472, sd=0.000046643)
c = rnorm(1000, mean=-0.000252, sd=0.00010862)
y = rnorm(1000, mean=16.53, sd=7.569)

x = (a+c*y)/-b

Where the means of a, b, and c are the parameter estimates from the regression and the SDs are their associated standard errors. The mean and SD of y are from summary statistics of the data.

What I am unsure of, is how to account for any covariance between any of a, b, and c? Or if there would be any issue to begin with?

P.S. I doubt the title of my question accurately reflects my question so it would be helpful if somebody could better characterize it

Answer 1

I see two viable options in your case (there is probably many more).

Option 1:

Since you have data, you could use bootstrap to resample a,b,c,x,y instead of assuming that they are independent with their respectives means and standard deviations. This will preserve the properties of your data (covariances as you want, but distributions also).

Option 2:

If you can not use bootstrap, you can also use the variance matrix of the estimator, which is $\text{var} (\hat{\beta})=\sigma^2 (X^\prime X)^{-1}$, where is $X$ is your predictors and $\sigma^2$ is the variance of residual error in your outcome. The diagonal elements are the same as you use as sd, but will include covariances also. You can then use this matrix to generate data as you have carried out. Using package MASS in R should easily do the job.