I would like to simulate a multiple linear regression model using R.
If I have the skewness and kurtosis for the residuals, how can I do that?
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It only takes a minute to sign up.Sign up to join this community
Although skewness and kurtosis do not determine the distribution of the residuals, you can choose a distribution (within mathematical constraints) that has zero mean, unit variance, and achieves the desired skewness and kurtosis (and then later rescale and shift it if desired). Perhaps the best known family of such distributions is the Pearson system. If such a distribution would be a suitable model for your residuals, then the simulation can proceed by:
Choosing values of all explanatory variables $x_1, \ldots, x_p$ and forming them by columns into a matrix $X$.
Choosing the coefficients $\beta_1, \ldots, \beta_p$.
Generating iid variates $\varepsilon$ from the Pearson distribution with the given skewness and kurtosis (and mean and variance).
Computing the dependent variable as $y = X\beta + \varepsilon$.
It is not guaranteed that the residuals from the model $y \sim X$ will have the given moments! The properties of the residuals depend on the method of fitting, among other things. But if the residuals are not terribly skew and are generated from a zero-mean distribution, then their distribution will generally cleave to the desired one.
Here is an
R implementation and illustration. It uses the PearsonDS library to generate the errors.
library(PearsonDS) set.seed(17) n <- 320 beta <- c(-1,2,3) moments <- c(0, 1, -1, 5) x <- cbind(rep(1, n), matrix(rnorm((length(beta)-1)*n), nrow=n)) e <- rpearson(n, moments=moments) y <- x %*% beta + e fit <- lm(y ~ x - 1) # The constant is included in `x` already par(mfrow=c(1,3)) hist(e, main="Errors") hist(residuals(fit), main="Residuals") plot(residuals(fit), e, main="Residuals vs. Errors")
The visual similarity of the histograms indicates that the randomly generated errors and the OLS residuals have qualitatively the same distribution (with a slightly negative skewness of $-1$ and slightly positive excess kurtosis of $5-3=2$). The tight scatterplot confirms this.