# 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. Mar 7, 2012 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, 2012 at 15:26
• I just add a real example, which should be helpful. thanks
– Tu.2
Mar 7, 2012 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)))


Bayesian networks (BNs) are commonly used in the context you describe. As a generative model, a BN would allow you to represent the statistical dependencies between your domain variables, which in your case can be subgrouped as 1) pre-treatment, 2) treatment, and 3) post-treatment variables. You can train your model on your existing patient data, and then enter evidence (fill in observed values) for a specific patient to investigate how the observed values affect other variables (including those you labelled as outcome, i.e. post-treatment.)

One neat trick is that you can actually asses the effect of different treatment types on your outcome variables. This is called an intervention. If interested, we have a relevant paper here.