data(DT_btdata)
DT <- DT_btdata
head(DT)
# \donttest{
############# sommer ################
if(requireNamespace("sommer")){
library(sommer)
mix4 <- mmes(tarsus ~ sex,
random = ~ dam + fosternest,
rcov=~units,
data = DT)
summary(mix4)$varcomp
# MULTI-TRAIT EXAMPLE
traits <- c("tarsus","back","hatchdate")
DT[,traits] <- apply(DT[,traits],2,scale)
DTL <- reshape(DT[,c("animal","sex","dam","fosternest", traits)],
idvar = c("animal","sex","dam","fosternest"),
varying = traits,
v.names = "value", direction = "long",
timevar = "trait", times = traits )
DTL <- DTL[with(DTL, order(trait,animal)), ]
head(DTL)
mix3 <- mmes(value ~ trait:sex - 1, henderson=TRUE,
random = ~ vsm(usm(trait),ism(dam)) +
vsm(usm(trait), ism(fosternest)),
rcov= ~ vsm(dsm(trait),ism(units)),
data = DTL)
summary(mix3)$varcomp
#### calculate the genetic correlation
cov2cor(mix3$theta[[1]])
cov2cor(mix3$theta[[2]])
}
############# lme4breeding ################
if(requireNamespace("lme4breeding")){
library(lme4breeding)
mix4 <- lmeb(tarsus ~ sex + (1|dam) + (1|fosternest),
data = DT)
vc <- VarCorr(mix4); print(vc,comp=c("Variance"))
sigma(mix4)^2 # error variance
BLUP <- ranef(mix4, condVar=TRUE)
PEV <- lapply(BLUP, function(x){attr(x, which="postVar")}) # take sqrt() for SEs
### multi-trait model
traits <- c("tarsus" ,"back", "hatchdate")
for(iTrait in traits){DT[,iTrait] <- scale(DT[,iTrait])}
DTL <- reshape(DT[,c("animal", traits)], idvar = "animal", varying = traits,
v.names = "value", direction = "long",
timevar = "trait", times = traits )
DTL <- merge(DTL, unique(DT[,c("animal","dam","fosternest","sex")]),
by="animal", all.x = TRUE)
DTL <- DTL[with(DTL, order(trait)), ]
head(DTL)
system.time(
mix <- lmeb(value ~ (0+trait|dam),
# relmat = list(geno=A),
# rotation = TRUE,
data=DTL)
)
vc <- VarCorr(mix); print(vc,comp=c("Variance"))
cov2cor(vc$dam)
sigma(mix)^2 # error variance
}
############# end ################
# }
Run the code above in your browser using DataLab