data(DT_rice)
DT <- DT_rice
GT <- GT_rice
GTn <- GTn_rice
head(DT)
M <- atcg1234(GT)
# \donttest{
##############################################
############### sommer. #################
##############################################
if(requireNamespace("sommer")){
library(sommer)
A <- A.mat(M$M)
mix <- mmes(Protein.content~1,
random = ~vsm(ism(geno), Gu=A) + geno,
rcov=~units,
data=DT)
summary(mix)$varcomp
# if using henderson=TRUE provide Gu as inverse
Ai <- solve(A + diag(1e-6,ncol(A),ncol(A)))
Ai <- as(as(as( Ai, "dMatrix"), "generalMatrix"), "CsparseMatrix")
attr(Ai, 'inverse')=TRUE
## MULTI-TRAIT MODEL
## reshape in long format the dataset
traits <- c("Protein.content","Flag.leaf.length")
DTL <- reshape(DT[,c("geno", traits)], idvar = "geno", varying = traits,
v.names = "value", direction = "long",
timevar = "trait", times = traits )
DTL <- DTL[with(DTL, order(trait)), ]
head(DTL)
M <- DTLM <- atcg1234(GT)
A <- A.mat(M$M)
mix <- mmes(value~trait,
random = ~vsm(usm(trait),ism(geno), Gu=A) ,
rcov=~vsm(dsm(trait), ism(units)),
data=DTL)
summary(mix)$varcomp
cov2cor(mix$theta$`vsm(usm(trait), ism(geno), Gu = A`)
}
##############################################
############### lme4breeding #################
##############################################
if(requireNamespace("lme4breeding")){
library(lme4breeding)
### univariate model
A <- A.matr(M$M)
A <- A + diag(1e-4, ncol(A), ncol(A))
mix <- lmeb(Protein.content ~ (1|geno),
relmat = list(geno=A),
data=DT)
vc <- VarCorr(mix); print(vc,comp=c("Variance"))
sigma(mix)^2 # error variance
### multi-trait model
traits <- c("Flowering.time.at.Arkansas" ,"Seed.volume", "Protein.content")
for(iTrait in traits){DT[,iTrait] <- scale(DT[,iTrait])}
DTL <- reshape(DT[,c("geno", traits)], idvar = "geno", varying = traits,
v.names = "value", direction = "long",
timevar = "trait", times = traits )
DTL <- DTL[with(DTL, order(trait, geno)), ]
head(DTL)
system.time(
mix <- lmeb(value ~ (0+trait|geno),
relmat = list(geno=A),
rotation = TRUE,
data=DTL)
)
vc <- VarCorr(mix); print(vc,comp=c("Variance"))
vc$geno
sigma(mix)^2 # error variance
BLUP <- ranef(mix, condVar=TRUE)
condVAR <- lapply(BLUP, function(x){attr(x, which="postVar")}) # take sqrt() for SEs
}
# }
Run the code above in your browser using DataLab