data(DT_cornhybrids)
DT <- DT_cornhybrids
GT <- GT_cornhybrids
A <- GT
K1 <- A[levels(DT$GCA1), levels(DT$GCA1)]; dim(K1)
K2 <- A[levels(DT$GCA2), levels(DT$GCA2)]; dim(K2)
# \donttest{
S <- kronecker(K1, K2) ; dim(S)
rownames(S) <- colnames(S) <- levels(DT$SCA)
# }
# \donttest{
################ sommer ###################
if(requireNamespace("sommer")){
library(sommer)
ans <- mmes(Yield ~ Location,
random = ~ vsm(ism(GCA1),Gu=K1) + vsm(ism(GCA2),Gu=K2), # + vsm(ism(SCA),Gu=S),
rcov=~units,
data=DT)
summary(ans)$varcomp
## mmec uses the inverse of the relationship matrix
K1i <- solve(K1 + diag(1e-4,ncol(K1),ncol(K1)))
K1i <- as(as(as( K1i, "dMatrix"), "generalMatrix"), "CsparseMatrix")
attr(K1i, 'inverse')=TRUE
K2i <- solve(K2 + diag(1e-4,ncol(K2),ncol(K2)))
K2i <- as(as(as( K2i, "dMatrix"), "generalMatrix"), "CsparseMatrix")
attr(K2i, 'inverse')=TRUE
Si <- solve(S + diag(1e-4,ncol(S),ncol(S)))
Si <- as(as(as( Si, "dMatrix"), "generalMatrix"), "CsparseMatrix")
attr(Si, 'inverse')=TRUE
ans2 <- mmes(Yield ~ Location,
random = ~ vsm(ism(GCA1),Gu=K1i) + vsm(ism(GCA2),Gu=K2i), # + vsm(ism(SCA),Gu=Si),
henderson=TRUE,
rcov=~units,
data=DT)
summary(ans2)$varcomp
}
################ lme4breeding ###################
if(requireNamespace("lme4breeding")){
library(lme4breeding)
ans <- lmeb(Yield ~ Location + (1| GCA1) + (1|GCA2),
relmat = list(GCA1=K1,
GCA2=K2),
data=DT)
vc <- VarCorr(ans); print(vc,comp=c("Variance"))
sigma(ans)^2 # error variance
BLUP <- ranef(ans, condVar=TRUE)
condVAR <- lapply(BLUP, function(x){attr(x, which="postVar")}) # take sqrt() for SEs
}
# }
Run the code above in your browser using DataLab