[R-sig-ME] Calculating variance-covariance matrices of random effects

[David A Springate]

Dear Everyone,

I am trying to calculate an additive genetic variance-covariance matrix
for a group of  t traits.
As an example, here is a toy dataset:

genotype = factor(rep(seq(101,125),each = 5))
trait1 = sapply(c(rep(seq(106,130),each = 5)), function(x) x + rnorm(1,
sd = 5))
trait2 = sapply(c(rep(seq(130,106),each = 5)), function(x) x + rnorm(1,
sd = 5))
trait3 = sapply(c(rep(runif(25,105,130),each = 5)), function(x) x +
rnorm(1, sd = 5))
trait4 = sapply(c(rep(seq(111,160,2),each = 5)), function(x) x +
rnorm(1, sd = 11))

Ultimately I want a t x t matrix, G, with among-genotype variance
components for all traits along the diagonal and among-genotype
covariance components for all traits on the off-diagonals. 

In this example, I would expect it to look something similar to the
output from:
traitmeans = sapply(list(trait1,trait2,trait3,trait4),
                    function(x) tapply(x,genotype,mean))
cov(traitmeans)

which, in my case is...
          [,1]        [,2]       [,3]        [,4]
[1,]  51.84631  -53.776108 -10.151024   92.786823
[2,] -53.77611   66.836473   8.126729 -108.332039
[3,] -10.15102    8.126729  63.895903   -1.148921
[4,]  92.78682 -108.332039  -1.148921  203.643692


To do this I think I need to do a random effects analysis of covariance
model for each trait and get the variances/covariances from @REmat.
I have tried:

mod1 = lmer(trait1 ~ (trait2 + trait3 + trait4|genotype))
mod2 = lmer(trait2 ~ (trait1 + trait3 + trait4|genotype))
mod3 = lmer(trait3 ~ (trait1 + trait2 + trait4|genotype))
mod4 = lmer(trait4 ~ (trait1 + trait2 + trait3|genotype))
 but the variances look completely wrong
e.g. calling...
summary(mod1) @ REmat   #gives

Groups     Name          Variance     Std.Dev.
Corr                      
 "genotype" "(Intercept)" "2.9504e+02" "17.176685" ""       ""
""      
 ""         "trait2"      "5.6112e-03" " 0.074908" "-1.000" ""
""      
 ""         "trait3"      "1.8227e-02" " 0.135006" "-0.998" " 0.998"
""      
 ""         "trait4"      "3.6162e-03" " 0.060135" " 0.577" "-0.577"
"-0.624"
 "Residual" ""            "2.8417e+01" " 5.330789" ""       ""
""      
I was expecting the first 4 rows in the variance column to be the first
column in my G matrix
The results from VarCorr(mod1) are no more enlightening, so I assume my
model must be wrong.
Any help would be much appreciated, as well as any suggestions of papers
using R for calculating variance components / G-matrices / multivariate
selection in quantitative genetics - they seem pretty thin on the
ground!

David