# How to simulate repeated measures multivariate outcomes in R?

@whuber has demonstrated how to simulate multivariate outcomes ($y_1$, $y_2$, and $y_3$) for one time point.

As we know, longitudinal data often occur in medical studies. My question is how to simulate repeated measures multivariate outcomes in R? For example, we repeatedly measure $y_1$, $y_2$, and $y_3$ at 5 various time points for two different treatment groups.

-

Use the rmvnorm() function, It takes 3 arguments: the variance covariance matrix, the means and the number of rows.

The sigma will have 3*5=15 rows and columns. One for each observation of each variable. There are many ways of setting these 15^2 parameters(ar, bilateral symmetry, unstructured...). However you fill in this matrix be aware of the assumptions, particularly when you set a correlation/covariance to zero, or when you set two variances to be equal. For a starting point a sigma matrix might might look something like this:

 sigma=matrix(c(
#y1             y2             y3
3 ,.5, 0, 0, 0, 0, 0, 0, 0, 0,.5,.2, 0, 0, 0,
.5, 3,.5, 0, 0, 0, 0, 0, 0, 0,.2,.5,.2, 0, 0,
0 ,.5, 3,.5, 0, 0, 0, 0, 0, 0, 0,.2,.5,.2, 0,
0 , 0,.5, 3,.5, 0, 0, 0, 0, 0, 0, 0,.2,.5,.2,
0 , 0, 0,.5, 3, 0, 0, 0, 0, 0, 0, 0, 0,.2,.5,
0 ,0 ,0 ,0 , 0, 3,.5, 0, 0, 0, 0, 0, 0, 0, 0,
0 ,0 ,0 ,0 ,0 ,.5, 3,.5, 0, 0, 0, 0, 0, 0, 0,
0 ,0 ,0 ,0 ,0 ,0 ,.5, 3,.5, 0, 0, 0, 0, 0, 0,
0 ,0 ,0 ,0 ,0 ,0 ,0 ,.5, 3,.5, 0, 0, 0, 0, 0,
0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,.5, 3, 0, 0, 0, 0, 0,
.5,.2,0 ,0 ,0 ,0 ,0 ,0 ,0 , 0, 3,.5, 0, 0, 0,
.2,.5,.2,0 ,0 ,0 ,0 ,0 ,0 ,0 ,.5, 3,.5, 0, 0,
0 ,.2,.5,.2,0 ,0 ,0 ,0 ,0 ,0 ,0 ,.5, 3,.5, 0,
0 ,0 ,.2,.5,.2,0 ,0 ,0 ,0 ,0 ,0 ,0 ,.5, 3,.5,
0 ,0 ,0 ,.2,.5,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,.5, 3

),15,15)


So the sigma[1,12] is .2 and that means that the covariance between the first observation of Y1 and the 2nd observation of Y3 is .2, conditional on all the other 13 variables. The diagonal rows do not all have to be the same number: that is a simplifying assumption that I made. Sometimes it makes sense, sometimes it doesn't. In general it means the correlation between a 3rd observation and a 4th is the same as the correlation between a 1st and a second.

You also need means. It could be as simple as

 meanTreat=c(1:5,51:55,101:105)
meanControl=c(1,1,1,1,1,50,50,50,50,50,100,100,100,100,100)


Here the first 5 are the means for the 5 observations of Y1, ... , the last 5 are the observations of Y3

then get 2000 observation of your data with:

sampleT=rmvnorm(1000,meanTreat,sigma)
sampleC=rmvnorm(1000,meanControl,sigma)
sample=data.frame(cbind(sampleT,sampleC) )
sample$group=c(rep("Treat",1000),rep("Control",1000) ) colnames(sample)=c("Y11","Y12","Y13","Y14","Y15", "Y21","Y22","Y23","Y24","Y25", "Y31","Y32","Y33","Y34","Y35")  Where Y11 is the 1st observation of Y1,...,Y15 is the 5th obs of Y1... - To create matrices as in the first example, Seth, try this: n <- 3*5; sigma <- diag(1, nrow=n, ncol=n); sigma[rbind(cbind(1:n-1,1:n),cbind(1:n,1:n-1))] <- 1/2. A similar approach will generate the second example. However, they have a common problem: you have lost the covariances among the$y$'s during each period--these matrices do not reflect a repeated measures structure. – whuber Apr 17 '12 at 17:29 @whuber your syntax is helpful, but is different from what I wrote. I think the difference is a bit important. I think of what I wrote as AR(1) and you have a entry in the cross correlation between the last observation of one variable and the first observation of the next variable. In other words I think sigma[5,6] should be 0. – Seth Apr 17 '12 at 17:55 Ah, now I see what you are doing: you are creating three AR(1) series in the first example. I missed this because I believe the OP is concerned also about correlations between the series: that's what is meant by "multivariate outcomes." From that point of view it seems like you want to view these correlation matrices as being$5$by$5$block matrices with each entry a$3$by$3$matrix, rather than as a$3$by$3$block matrix with$5$by$5$blocks. – whuber Apr 17 '12 at 18:58 I thought my second sigma was a simple example of allowing the variance between Y1 and Y3 to be positive. I edited the answer a bit to make clearer that the matrix is there to be configured depending on the data generating process. There are definitely many ways to skin this cat. – Seth Apr 17 '12 at 19:06 Fair enough, but your approach creates difficulties, because it's not trivial to combine the multivariate correlation among the$y_i\$ with an AR model. E.g., were you aware that the second matrix fails to be positive definite? (The "variance" of c(-102, 177, -204, 177, -102, 0, 0, 0, 0, 0, 102, -177, 204, -177, 102) is negative.) – whuber Apr 17 '12 at 19:11

To generate multivariate normal data with a specified correlation structure, you need to construct the variance covariance matrix and calculate its Cholesky decomposition using the chol function. The product of the Cholesky decomposition of the desired vcov matrix and independent random normal vectors of observations will yield random normal data with that variance covariance matrix.

v <- matrix(c(2,.3,.3,2), 2)
cv <- chol(v)

o <- replicate(1000, {
y <- cv %*% matrix(rnorm(100),2)

v1 <- var(y[1,])
v2 <- var(y[2,])
v3 <- cov(y[1,], y[2,])

return(c(v1,v2,v3))
})

## MCMC means should estimate components of v
rowMeans(o)

-