title: Simulating data from multivariate distribution in R based on Winbugs/JAGS script
body: <p>I am trying to simulate data, based on part of a JAGS/Winbugs script. The script comes from Eaves &amp; Erkanli (2003, see, <a href="http://psych.colorado.edu/~carey/pdffiles/mcmc_eaves.pdf" rel="nofollow">http://psych.colorado.edu/~carey/pdffiles/mcmc_eaves.pdf</a>, page 295-296). </p>
  
  <p>The (part of) the script I want to base my simulations on is as follows (different variable names than in the original paper): </p>
  
  <pre><code> for(fam in 1 : nmz ){
      a2mz[fam, 1:N] ~ dmnorm(mu[1:N], tau.a[1:N, 1:N])
      a1mz[fam, 1:N] ~ dmnorm(a2mz[fam, 1:N], tau.a[1:N, 1:N])
   }
  
   #Prior
   tau.a[1:N, 1:N] ~ dwish(omega.g[,], N) 
  </code></pre>
  
  <p>I want to simulate data in R for the parameters a2mz and a1mz as given in the script above.</p>
  
  <p>So basically, I want to simualte data from -N- (e.g. =  3) multivariate distributions with -fam- (e.g. 10) persons with sigma tau.a. </p>
  
  <p>To make this more illustrative: The purpose is to simulate genetic effects for -fam- (e.g. 10) families. The genetic effect is the same for each family (e.g. monozygotic twins), with a variance of tau.a (e.g. 0.5). Of these genetic effects, 3 'versions' (3 multivariate distributions) have to be simulated. </p>
  
  <p>What I tried in R to simulate the data as given in the JAGS/Winbugs script is as follows: </p>
  
  <pre><code> library(MASS)
   nmz = 10 #number of families, here e.g. 10
   var_a = 0.5 #tau.g in the script
  
   a2_mz &lt;- mvrnorm(3, mu = rep(0, nmz), Sigma = diag(nmz)*var_a)
  </code></pre>
  
  <p>This simulates data for the a2mz parameter as referred to in the JAGS/Winbugs script above:</p>
  
  <pre><code> &gt; print(t(a2_mz))
          [,1]       [,2]        [,3]
    [1,] -1.1563683 -0.4478091 -0.15037563
    [2,]  0.5673873 -0.7052487  0.44377336
    [3,]  0.2560446  0.9901964 -0.65463341
    [4,] -0.8366952  0.4924839 -0.56891991
    [5,]  0.7343780  0.5429955  0.87529201
    [6,]  0.5592868 -0.3899988 -0.33709105
    [7,] -1.8233663 -0.7149141 -0.18153049
    [8,] -0.8213804 -1.4397075 -0.09159725
    [9,] -0.7002797 -0.3996970 -0.29142215
    [10,]  1.1084067  0.3884869 -0.46207940 
  </code></pre>
  
  <p>However, when I then try to use these data to simulate data for the a1mz (third line of the JAGS/Winbugs) script, then something goes wrong and I am not sure what:</p>
  
  <pre><code> a1_mz &lt;- mvrnorm(3, mu = t(a2_mz), Sigma = c(diag(nmz)*var_a, diag(nmz)*var_a,     diag(nmz)*var_a))
  </code></pre>
  
  <p>This results in the error: 
  Error in eigen(Sigma, symmetric = TRUE, EISPACK = EISPACK) : 
    non-square matrix in 'eigen'</p>
  
  <p>Can anyone give me any hints or tips on what I am doing wrong? </p>
  
  <p>Many thanks, 
  Best regards, 
  inga </p>