Simulating non-normal correlated data for Bayesian regression

Question

I'm interested in generating data for three separate datasets where each contains three IVs and a single DV that are correlated with one another based on meta-analytic data. For example, I would like the correlation matrix between the variables to look as follows based on the meta-analytic data:

IV1 IV2 IV3 DV 

IV1 1
IV2 0.31 1
IV3 0.14 0.21 1
DV 0.03 0.18 0.18 1

Another nuanced bit to the simulated data is that I would like those relationships to hold across the datasets, but I want to manipulate the distribution of the dependent variable. In the first dataset, the distribution of the DV would be normal, the second dataset the DV would have a moderate, positive skew, and the third dataset the DV would follow a power-law distribution.

I haven't simulated data before so I am very naive to the process, so any help would be greatly appreciated.

Thanks in advance!

Answer 1

The challenge is creating the intended correlations while, at the same time, achieving the intended non-normal marginals that you mention. Here are two options with references and links:

1. Headrick (2002)

This approach is very fast, but it cannot create all possible distribution shapes. Its usefulness for you will depend on how extreme your power law distribution is (if $\alpha$ is too low, I doubt this method will work). You can find the details in the references below, and I'll summarize it briefly. The basic idea is to simulate standard normal variables, rotate them using intermediate correlations, and then construct your intended variables each as a polynomial of these z-scores. For example, $Y=c_0+c_1z+c_2z^2+c_3z^3+c_4z^4+c_5z^5$. The most difficult part is finding the constants, which is typically done in Mathematica. A tutorial can be found in Headrick et al. (2007).

Headrick, T. C. (2002). Fast fifth-order polynomial transforms for generating univariate and multivariate nonnormal distributions. Computational Statistics & Data Analysis, 40(4), 685-711.

Headrick, T. C., Sheng, Y., & Hodis, F. A. (2007). Numerical Computing and Graphics for the Power Method Transformation Using Mathematica. Journal of Statistical Software, 19(3).

2. Ruscio & Kaczetow (2008)

This is a slower, iterative algorithm. Their article is behind a paywall, but their R code is freely available here. I'll summarize the bivariate case for the sake of simplicity, though it can work with multivariate problems, too. Uncorrelated $X_o$ and $Y_o$ are generated with any shapes (e.g. bimodal). Then, $X_1$ and $Y_1$ are generated as bivariate normal with an intermediate correlation. $X_1$ and $Y_1$ are replaced by $X_0$ and $Y_0$ in a rank-preserving fashion. Adjust the intermediate correlation higher or lower depending on the observed correlation between $X_1$ and $Y_1$. Rinse and repeat until the observed correlation is sufficiently close to the target correlation.

Ruscio, J., & Kaczetow, W. (2008). Simulating multivariate nonnormal data using an iterative algorithm. Multivariate Behavioral Research, 43, 355–381. doi:10.1080/00273170802285693