[R-sig-ME] Calculating variance-covariance matrices of random effects
David A Springate
david.springate at postgrad.manchester.ac.uk
Thu Feb 10 14:51:26 CET 2011
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
More information about the R-sig-mixed-models
mailing list