# NOT RUN {
####=========================================####
#### For CRAN time limitations most lines in the
#### examples are silenced with one '#' mark,
#### remove them and run the examples using
#### command + shift + C |OR| control + shift + C
####=========================================####
####=========================================####
####=========================================####
#### EXAMPLE 1
#### breeding values with 3 variance components
#### for hybrid prediction
####=========================================####
####=========================================####
data(cornHybrid)
hybrid2 <- cornHybrid$hybrid # extract cross data
A <- cornHybrid$K
y <- hybrid2$Yield
####=========================================####
#### Realized IBS relationships for set of parents 1
####=========================================####
K1 <- A[levels(hybrid2$GCA1), levels(hybrid2$GCA1)]; dim(K1)
####=========================================####
#### Realized IBS relationships for set of parents 2
####=========================================####
K2 <- A[levels(hybrid2$GCA2), levels(hybrid2$GCA2)]; dim(K2)
####=========================================####
#### SCA relationship matrix (kronecker)
####=========================================####
S <- kronecker(K1, K2, make.dimnames=TRUE) ; dim(S)
#head(hybrid2)
#ans <- mmer2(Yield ~ Location, random = ~ g(GCA1) + g(GCA2) + g(SCA),
# G=list(GCA1=K1, GCA2=K2, SCA=S),data=hybrid2)
#ans$var.comp
#summary(ans)
#### you can specify heterogeneus variances
# ans <- mmer2(fixed=Yield ~ 1, random = ~ GCA2 + at(Location):GCA2,
# data=hybrid2)
# summary(ans)
########################################################
########################################################
########################################################
########################################################
########################################################
########################################################
####=========================================####
####=========================================####
#### EXAMPLE 2
#### breeding values for 1 variance component
#### with variance covariance structure
####=========================================####
####=========================================####
data(CPdata)
CPpheno <- CPdata$pheno
CPgeno <- CPdata$geno
#### create the variance-covariance matrix
A <- A.mat(CPgeno)
#### look at the data and fit the model
head(CPpheno)
mix1 <- mmer2(Yield~1,random=~g(id), G=list(id=A), data=CPpheno)
summary(mix1)
########################################################
########################################################
########################################################
########################################################
########################################################
########################################################
####=========================================####
#### EXAMPLE 3
#### comparison with lmer, install 'lme4'
#### and run the code below
####=========================================####
#### lmer cannot use var-cov matrices so we will not
#### use them in this comparison example
#library(lme4)
#library(sommer)
#data(cornHybrid)
#hybrid2 <- cornHybrid$hybrid
#fm1 <- lmer(Yield ~ Location + (1|GCA1) + (1|GCA2) + (1|SCA),
# data=hybrid2 )
#out <- mmer2(Yield ~ Location, random = ~ GCA1 + GCA2 + SCA,
# data=hybrid2)
#summary(fm1)
#summary(out)
#### same BLUPs for GCA1, GCA2, SCA than lme4
#plot(out$cond.residuals, residuals(fm1))
#plot(out$u.hat[[1]], ranef(fm1)$GCA1[,1])
#plot(out$u.hat[[2]], ranef(fm1)$GCA2[,1])
#vv=which(abs(out$u.hat[[3]]) > 0)
#plot(out$u.hat[[3]][vv,], ranef(fm1)$SCA[,1])
#### a more complex model specifying which locations
#out2 <- mmer2(Yield ~ 1, random = ~ at(Location,c("3","4")):GCA2
# + at(Location,c("3","4")):SCA,
# data=hybrid2)
#summary(out2)
########################################################
########################################################
########################################################
########################################################
########################################################
########################################################
####=========================================####
#### EXAMPLE 4
#### comparison with lmer, install 'lme4'
#### and run the code below
####=========================================####
#library(lme4)
#data(sleepstudy)
#fm1 <- lmer(Reaction ~ Days + (1 | Subject), sleepstudy)
#summary(fm1)
#out <- mmer2(Reaction ~ Days, random = ~ Subject, data=sleepstudy)
#summary(out)
# plot(out$cond.residuals, residuals(fm1))
# plot(out$u.hat[[1]], ranef(fm1)[[1]][,1])
### specific variances for Subjects at Days
# sleepstudy$Days <- as.factor(sleepstudy$Days)
# out <- mmer2(Reaction ~ 1, random = ~ Subject + at(Days):Subject,
# data=sleepstudy)
# summary(out)
########################################################
########################################################
########################################################
########################################################
########################################################
########################################################
####=========================================####
#### EXAMPLE 5
#### Multivariate model example
####=========================================####
# data(FDdata)
# head(FDdata)
#
# mix <- mmer2(cbind(stems,pods,seeds)~1,
# random=~female+male, MVM=TRUE,data=FDdata)
# summary(mix)
# #### genetic variance covariance
# gvc <- mix$var.comp$female
# #### extract variances (diagonals) and get standard deviations
# sd.gvc <- as.matrix(sqrt(diag(gvc)))
# #### get possible products sd(Vgi) * sd(Vgi')
# prod.sd <- sd.gvc %*% t(sd.gvc)
# #### genetic correlations cov(gi,gi')/[sd(Vgi) * sd(Vgi')]
# (gen.cor <- gvc/prod.sd)
# #### pods and seeds have a strong negative genetic covariance (-.79)
# #### more pods, less seeds
########################################################
########################################################
########################################################
########################################################
########################################################
########################################################
####=========================================####
#### EXAMPLE 6
#### More on modeling
####=========================================####
# library(lme4)
# data(sleepstudy)
# head(sleepstudy)
# # try lme4
# (fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy))
# summary(fm1)
# # try sommer
# out <- mmer2(Reaction ~ Days, random = ~ Subject + at(Days):Subject, data=sleepstudy)
# summary(out)
# plot(out$cond.residuals, residuals(fm1))
# plot(out$u.hat$Subject, ranef(fm1)[[1]][,1])
# # now consider Days as factor by creating a new indicator variable
# sleepstudy$Days2 <- as.factor(sleepstudy$Days)
# out2 <- mmer2(Reaction ~ Days, random = ~ Subject + at(Days2):Subject, data=sleepstudy)
# summary(out2)
# plot(out2$cond.residuals, residuals(fm1))
# plot(out2$u.hat$Subject, ranef(fm1)[[1]][,1])
# }
Run the code above in your browser using DataLab