# Problem when creating matrix of values based on covariance matrix

I want to simulate a data set with similar covariance structure as my observed data (which is a SNP by gene p-value matrix, dim ~600k*8368), and have calculated a covariance matrix (dimensions 8368*8368). So far I've tried two approaches to simulating data (all in R):

rmvnorm from the mvtnorm package, and by using the Cholesky decomposition.

The problem is that the resulting covariance structure in my simulated data is very different from the original supplied covariance vector (but the code used is accurate as far as I can tell). Lets just look at some of the values:

    > cov8[1:4,1:4] # covariance of simulated data
X1          X2         X3         X4
X1 34515296.00    99956.69   369538.1  1749086.6
X2    99956.69 34515296.00  2145289.9  -624961.1
X3   369538.08  2145289.93 34515296.0  -163716.5
X4  1749086.62  -624961.09  -163716.5 34515296.0
> CEUcovar[1:4,1:4]
[,1]         [,2]          [,3]         [,4]
[1,] 0.1873402987  0.001837229  0.0009009272  0.010324521
[2,] 0.0018372286  0.188665853  0.0124216535 -0.001755035
[3,] 0.0009009272  0.012421654  0.1867835412 -0.000142395
[4,] 0.0103245214 -0.001755035 -0.0001423950  0.192883488


So the distribution of the observed covariance is very narrow compared to the simulated data.

None of the eigenvalues of the observed covariance matrix are negative, and it appears to be a positive definite matrix. Here is what I did to create the simulated data (very straightforward):

    Chol <- chol(CEUcovar)
Z <- matrix(rnorm(20351 * 8368), 8368)
X <- t(Chol) %*% Z
sample8 <- data.frame(as.matrix(t(X)))
> dim(sample8)
[1] 20351  8368
cov8=cov(sample8,method='spearman')


(earlier I've also tried: sample8 <- rmvnorm(1000,mean=rep(0,ncol(CEUcovar)), sigma=CEUcovar, method="eigen") with as 'bad' results, much larger covariance values in the simulated data)

• Any ideas of WHY the simulated data have such a different covariance?
• Any experience with similar issues? Would be happy to supply the covariance matrix if anyone wants to give it a try.
• Is this an impossible task statistically?
• Any suggestions? Anything apparent that I left our or neglected?

I am aware that my supplied covariance matrix might differ from the true underlying covariance of the population. But my goal is to simulate a random data that captures (as much as possible) of the covariance that this observed data has, no matter how weird that covariance is.

Clarification - my supplied covariance vector has values in [-10^-2, 0.3] while the resulting data has covariance of [-800000,5000000]. This is not what I am expecting. I wonder why this is the case, and how this might be explained.

Any advice would be highly appreciated.

-
Welcome to the site, @Boo. Do you see this question as having statistical content, or only about how to get this done in R (ie a coding question)? If the latter, it really belongs on Stack Overflow, where it would also be more likely to get a satisfactory answer (although that doesn't make it a bad question). Please don't cross-post (SE strongly discourages this). If you think it should go to SO, just say so, & after a bit the moderators can migrate it for you. –  gung Aug 11 '12 at 15:11
Thanks gung. I'm believe I am using the correct code in R. So this might either be a statistical question, ie whether I am trying to do something that is simply not possible. Or the problem might be due to an issue that R has with this covariance matrix (and the question is posted on R-help stat.ethz.ch/pipermail/r-help/2012-August/321409.html). –  Boo Aug 11 '12 at 15:20
@Boo This might be related to the sample size and the speed of convergence of the estimation. You are trying to estimate p entries using $\approx 2.5p$ observations. –  user10525 Aug 11 '12 at 18:10
Thanks @Procrastinator. Do you mean because the original data and the simulated data are of diff sizes (1/30 of the size)? But I am simply using the covariance. –  Boo Aug 11 '12 at 18:12
@Boo (It is possible that I am not properly understanding the problem). Yes, I tried on lower dimensions using for example a diagonal covariance matrix of size $4\times 4$ and estimating it using $10-30$ observations. The estimation is in many cases misleading. In your case you are estimating $8368\times 8367/2$ entries using $20351\times 8368$ observations. –  user10525 Aug 11 '12 at 18:18