Simulate non-normal data for multiple regression In a working paper, Luo (2011) analyzes generating non-normal data, using Fleishman power method. But this approach just can be used for specific skewness and kurtosis. I would like to simulate non-normal data for larger values of skewness and kurtosis (for example, skewness = 3.5 and kurtosis = -4.0). My questions are: 1) how can I do this; 2) are there any R package that can generate such data; 3) how can I generate non-normal data for residuals of the regression model?
References
Luo, H. (2011). Generation of non-normal data – A study of Fleishman’s power method. [Working paper]. Uppsala University. Retrieved from www.diva-portal.org/smash/get/diva2:407995/FULLTEXT01.pdf
 A: How to generate non-normal data with specific skewness and kurtosis values?
The Fleishman power method (a.k.a. 3rd order power polynomial) is able to handle a wide range of skewness and kurtosis values, but not all of them.  Keep in mind that kurtosis is bounded below by skewness (see Relationship between skew and kurtosis in a sample).  The Fleishman method has trouble simulating distributions that are even close to this boundary.  A slightly more flexible method is Headrick and Sawilowsky's (1999) 5th order power polynomial method.  You can find extensive details about that method and related ones in Headrick's (2009) book.  These polynomial methods can work well for many smooth distributions, but they have trouble simulating uniform and most multimodal distributions.  For such distributions, you might have to generate some known family (e.g., just generate a uniform with runif()). 
Are there any R packages that can generate such data?
Generating the data should be easy using rnorm() and the equations in the cited articles.  The hard part is finding the constants to use for specific skewness and kurtosis values.  For the Fleishman 3rd order polynomial, Zopluoglu (2011) has R code.  For the 5th order power polynomial method, I'm not aware of a reliable solution in R, and I've used Mathematica for that.  You can find the Mathematica code here.
How can I generate non-normal data for residuals of the regression model?
Just add your non-normal data to the regression equation.  For example, let $E$ be your vector of non-normal residuals:
$Y_i = B_0 + B_1*X_i + E_i$
References:
Headrick, T. C. (2009). Statistical simulation: power method polynomials and other transformations. CRC Press.
Headrick, T. C., & Sawilowsky, S. S. (1999). Simulating correlated multivariate nonnormal distributions: Extending the Fleishman power method. Psychometrika, 64(1), 25-35.
Ruscio, J., & Kaczetow, W. (2008). Simulating multivariate nonnormal data using an iterative algorithm. Multivariate Behavioral Research, 43(3), 355-381.
Zopluoglu, C. (2011). Applications in R: Generating multivariate non-normal variables. http://www.tc.umn.edu/~zoplu001/resim/gennonnormal.pdf.
