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

Answer 2

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.