data(DT_example)
DT <- DT_example
A <- A_example
head(DT)
# \donttest{
################## sommer #####################
if(requireNamespace("sommer")){
library(sommer)
## Compound simmetry (CS) model
ans1 <- mmes(Yield~Env,
random= ~ Name + Env:Name,
rcov= ~ units,
data=DT)
summary(ans1)$varcomp
####===========================================####
#### Univariate heterogeneous variance models ####
####===========================================####
## Compound simmetry (CS) + Diagonal (DIAG) model
DT=DT[with(DT, order(Env)), ]
ans2 <- mmes(Yield~Env,
random= ~Name + vsm(dsm(Env),ism(Name)),
rcov= ~ vsm(dsm(Env),ism(units)),
data=DT)
summary(ans2)
####===========================================####
#### Multi-trait variance models ####
####===========================================####
# stack traits
traits <- c("Yield","Weight")
DT[,traits] <- apply(DT[,traits],2,scale)
DTL <- reshape(DT[,c("Name","Env","Block", traits)],
idvar = c("Name","Env","Block"),
varying = traits,
v.names = "value", direction = "long",
timevar = "trait", times = traits )
DTL <- DTL[with(DTL, order(trait,Env)), ]
head(DTL)
## model
ans1 <- mmes(value~ trait,
random= ~ vsm(usm(trait), ism(Name)),
rcov= ~ vsm(dsm(trait), ism(units)),
data=DTL)
summary(ans1)$varcomp
cov2cor(ans1$theta$`vsm(usm(trait), ism(Name`)
}
################## lme4breeding #####################
if(requireNamespace("lme4breeding")){
library(lme4breeding)
## Compound simmetry (CS) model
ans1 <- lmeb(Yield~Env + (1|Name) + (1|Env:Name),
data=DT)
vc <- VarCorr(ans1); print(vc,comp=c("Variance"))
BLUP <- ranef(ans1, condVar=TRUE)
condVAR <- lapply(BLUP, function(x){attr(x, which="postVar")}) # take sqrt() for SEs
## Main (M) + Diagonal (DIAG) model
## with relationship matrix
ansCSDG <- lmeb(Yield ~ Env + (Env||Name),
relmat = list(Name = A ),
data=DT)
vc <- VarCorr(ansCSDG); print(vc,comp=c("Variance"))
## Main (M) + Diagonal (DIAG) model
## with diagonal residuals
## with relationship matrix
ansCSDG <- lmeb(Yield ~ Env + (Env||Name) + (0+Env||unitsR),
relmat = list(Name = A ),
data=DT)
vc <- VarCorr(ansCSDG); print(vc,comp=c("Variance"))
sigma(ansCSDG)^2 # error variance
}
# }
Run the code above in your browser using DataLab