How to simulate multivariate outcomes in R?

Question

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.

Answer 1

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)))