8
$\begingroup$

For regression analysis, it is often useful to know the data generating process to check how the used method works. While it is fairly simple to do this for a simple linear regression, this is not the case when the dependent variable has to follow a specific distribution.

Consider a simple linear regression:

N    <- 100
x    <- rnorm(N)
beta <- 3 + 0.4*rnorm(N)
y    <- 1 + x * beta + .75*rnorm(N)

Is there any way to use the same approach but to get y be other than normal, say left skewed?

$\endgroup$
10
$\begingroup$

If I understand your question correctly, this is quite easy. You just need to decide what distribution you want your errors to have, and use the corresponding random generation function.

There are a number of skewed distributions, so you need to figure out which one you like. In addition, most skewed distributions (e.g., log normal, chi-squared, Gamma, Weibull, etc.) are right skewed, so some minor adaptations would be necessary (e.g., multiply by $-1$).

Here is an example modifying your code:

set.seed(5840)  # this makes the example exactly reproducible
N      <- 100
x      <- rnorm(N)
beta   <- 0.4
errors <- rlnorm(N, meanlog=0, sdlog=1)
errors <- -1*errors   # this makes them left skewed
errors <- errors - 1  # this centers the error distribution on 0
y      <- 1 + x*beta + errors

enter image description here

I should note at this point that regression does not make any assumptions about the distributions of $X$ or $Y$, only about the errors, $\varepsilon$ (see here: What if the residuals are normally distributed, but y is not?). Thus, that was the focus of my answer above.


Update: Here is a right-skewed version with the errors distributed as Weibull:

set.seed(5840)  # this makes the example exactly reproducible
N      <- 100
x      <- rnorm(N)
beta   <- 0.4
errors <- rweibull(N, shape=1.5, scale=1)
# errors <- -1*errors   # this makes them left skewed
errors <- errors - factorial(1/1.5)  # this centers the error distribution on 0
y      <- 1 + x*beta + errors

Weibull data are right skewed already, so we don't need to switch their direction (i.e., we drop the -1*errors part). Also, from the Wikipedia page for the Weibull distribution, we see that the mean of a Weibull should be: $E[W] = (1/{\rm shape})!$. We want to subtract that value from each of the errors so that the resulting error distribution is centered on $0$. That allows the structural part (i.e., 1 + x*beta) of your code to accurately reflect the structural part of the data generating process.

enter image description here

The ExGaussian distribution is the sum of a normal and an exponential. There is a function ?rexGAUS in the gamlss.dist package to generate these. I don't have that package, but you should be able to adapt my code above without too much difficulty. You could also generate a random normal variable (via rnorm()) and an exponential (via rexp()) and sum them quite easily. Just remember to subtract the population mean, $\mu + 1/\lambda$, from each error prior to adding the errors to the structural part of the data generating process. (Be careful not to subtract the sample mean, mean(errors), though!)


Some final, unrelated comments: Your example code in the question is somewhat muddled (meaning no offense). Because rnorm(N) generates data with mean=0 and sd=1 by default, 0.4*rnorm(N) will generate rnorm(N, mean=0, sd=0.4). Your code (and possibly your thinking) will be much clearer if you use the latter formulation. In addition, your code for beta seems confused. We generally think of the $\beta$ in a regression-type model as a parameter, not a random variable. That is, it is an unknown constant that governs the behavior of the data generating process, but the stochastic nature of the process is encapsulated by the errors. This isn't the way we think about it when we are working with multilevel models, and your code seems halfway between a standard regression model and the code for a multilevel regression model. Specifying your betas separately is a good idea for maintaining the conceptual clarity of the code, but for a standard regression model, you would just assign a single number to each beta (e.g., beta0 <- 1; beta1 <- .04).

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ thanks for your answer, to be more specific: I need to simulate response time data. Response times are typically right skewed and are approximately distributed like an Ex-Gaussian or Weibull. $\endgroup$ – beginneR May 24 '14 at 14:44
  • $\begingroup$ I'll work up such an example, but it will be the errors that are skewed, the marginal distribution of Y is incidental. Also, response time data should be right skewed (a few long RTs, but a solid floor), not left skewed (towards small values). $\endgroup$ – gung - Reinstate Monica May 24 '14 at 14:51
  • $\begingroup$ Yes right-skewed, I accidentally wrote "left-skewed" in my question. $\endgroup$ – beginneR May 24 '14 at 14:54
  • 1
    $\begingroup$ thanks for your extensive answers. Actually, I am simulating data for a multilevel model. I just used the code and adopted it for a simple linear regression scenario. $\endgroup$ – beginneR May 24 '14 at 15:53
  • $\begingroup$ one more question: do you know how to simulate repeated measures data using a regression model like the one above so that the repeated measures have a specific correlation structure (to account for sphericity)? Say there are N subjects each providing 20 trials in each of three treatment conditions where the correlation between the treatment conditions should be known in the data generating process. $\endgroup$ – beginneR May 25 '14 at 10:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.