How to simulate multivariate outcomes in R?

Most of situations, we only deal with one outcome/response variable such as $y = a + bx +\epsilon$. However, in some scenarios, especially in the clinical data, the outcome variables can be high-dimensional/multivariate. Such as $\mathsf{Y} = \beta{x} + \mathsf{\epsilon}$, where $\mathsf{Y}$ contains $Y_1$, $Y_2$ and $Y_3$ variables and these outcomes are all correlated. If $x$ represents receiving treatment (yes/no), how can I simulated this type of data in R?

A real life example, each patient receives one of 2 types of bypass surgeries and researchers measure each patients on pain, swelling, fatigue ...etc after the bypass surgery (each symptom rates from 0 to 10). I "assume" outcomes(symptom severities) are multivariate normal. Hope this real example can clarify my question. Many thanks in advance.

-
What distribution does ${\bf Y}$ have? If it's multivariate normal, have a look at the mnormt library in R. –  Macro Mar 7 '12 at 15:26
This question is rather broad because "multivariate data" covers a lot of ground. What specific application do you have in mind? –  whuber Mar 7 '12 at 15:26
I just add a real example, which should be helpful. thanks –  Tu.2 Mar 7 '12 at 16:19

Simulate multivariate normal values with mvtnorm::rmvnorm. It doesn't seem to work quite like the univariate random number generators, which allow you to specify vectors of parameters, but this limitation is straightforward to work around.

For example, consider the model

$$E(y_1,y_2,y_3) = (-1+x, 2x, 1-3x)$$

where $\mathbf{y}$ has a multivariate normal distribution and $\text{Var}(y_i)=1$, $\text{Cov}(y_1, y_2) = \text{Cov}(y_2, y_3) = 0.5$, and $\text{Cov}(y_1,y_3)=0$. Let's specify this covariance matrix in R:

sigma <- matrix(c(1,   0.5, 0,
0.5, 1,   0.5,
0,   0.5, 1  ), 3, 3)


To experiment, let's generate some data for this model by letting $x$ vary from $1$ through $10$, with three replications each time. We have to include constant terms, too:

data <- cbind(rep(1,10*3), rep(1:10,3))


The model determines the means:

beta <- matrix(c(-1,1,  0,2,  1,-3), 2, 3)
means <- data %*% beta


The workaround for generating multiple multivariate results is to use apply:

library(mvtnorm) # Contains rmvnorm
sample <- t(apply(means, 1, function(m) rmvnorm(1, mean=m, sigma=sigma)))

-
thanks. This is really helpful –  Tu.2 Mar 7 '12 at 19:52