sommer (version 3.1)

mmer2: mixed model equations in R

Description

Fits a multivariate or univariate linear mixed model by likelihood methods (REML). It has been implemented to work with a data frame in a formula-based fashion which will create the incidence and variance covariance matrices for the user and call the mmer function in the background. The optimization methods are; Direct-Inversion Newton-Raphson (NR) (Tunnicliffe 1989; Lee et al. 2015; Maier et al. 2015), and Efficient Mixed Model Association (EMMA) (Kang et al. 2008). These algorithms are intended to be used for problems of the type p > n. For problems with sparse covariance structures, or problems of the type n > p, the MME-based algorithms are faster and we recommend to shift to the use of MME-based software (i.e. lme4, breedR, or asreml-R).

The package also provides kernels to estimate additive (A.mat), dominance (D.mat), and epistatic (E.mat) relationship matrices to model covariances among genotypes for plant and animal breeding problems.

For a short tutorial of how to perform different genetic analysis with sommer please type:

vignette("sommer.start") or vignette("sommer")

Usage

mmer2(fixed, random, rcov, data, G=NULL, grouping=NULL, method="NR",
      init=NULL,iters=20,tolpar = 1e-06, tolparinv = 1e-06, 
      draw=FALSE, silent=FALSE, constraint=TRUE, EIGEND=FALSE,
      forced=NULL, complete=TRUE, restrained=NULL,
      na.method.X="exclude", na.method.Y="exclude", REML=TRUE)

Arguments

fixed

a formula specifying the response variable(s) and fixed effects, i.e:

Yield ~ Location for univariate models

cbind(Yield,color) ~ Location for multivariate models.

random

a formula specifying the name of the random effects, i.e. random= ~ genotype + year.

Useful functions can be used to fit heterogeneous variances and other special models:

at(.) can be used to specify heterogeneous variance, i.e. random=~at(Location,c("A","B")):ID

diag(.) can be used to specify a diagonal variance structure, i.e. random=~diag(Location):ID

us(.) can be used to specify a unstructured variance, i.e. random=~us(Location):ID

and(.) can be used to specify overlay of matrices, i.e. random=~male + and(female)

g(.) can be used to specify that a covariance structure for the random effect (not inverse) should be included, i.e. random=~g(male), G=list(male=Am), etc.

grp(.) can be used to specify customized incidence matrices for random effects, i.e. random=~male + grp(re), grouping=list(re=Zre)

For multivariate models the functions; us(trait) and diag(trait) can be used before the random effects to indicate the covariance structure among traits, i.e.:

us(trait):Geno or diag(trait):Geno

us(trait):at(Location):Geno or diag(trait):at(Location):Geno

us(trait):us(Location):Geno or diag(trait):us(Location):Geno

rcov

a formula specifying the name of the error term, i.e. rcov= ~ units or rcov=~at(.):units.

The functions that can be used to fit heterogeneous residual variances are:

at(.) can be used to specify heterogeneous variance, i.e. rcov=~at(Location):units

For multivariate models the functions; us(trait) and diag(trait) can be used before the residual term to indicate the covariance structure among traits, i.e.:

us(trait):units or diag(trait):units

us(trait):at(Location):units or diag(trait):at(Location):units

data

a data frame containing the variables specified in the formulas for response, fixed, and random effects.

G

a list containing the variance-covariance matrices for the random effects (not inverses), i.e. if the user specifies:

random=~ g(genotype) + g(year)

then the G argument should be used as:

G=list(genotype=M1, year=M2)

where M1 and M2 are the variance-covariance matrices for random effects 'genotype' and 'year' respectively. Opposite to asreml you provide the original var-covar matrix, not its inverse.

grouping

argument to provide customized random effects. This has to be provided as a list where each element of the list is an incidence matrix for a random effect specified in the random formula using the `grp()` function. For example a customized random effect 're' would specified as:

random=~ grp(re), grouping=list(re=Zre)

where Zre is the customized incidence matrix. Additional functions such as g() can be used and the variance covariance matrix is provided in the G argument as any other random effect.

method

this refers to the method or algorithm to be used for estimating variance components. Sommer now focus only in 2 algorithms; Direct-inversion Newton-Raphson NR (Tunnicliffe 1989; Gilmour et al. 1995; Lee et al. 2015), and EMMA efficient mixed model association (Kang et al. 2008).

init

initial values for the variance components. By default this is NULL and variance components are estimated by the method selected, but in case the user want to provide initial values this argument is functional. It has to be provided as a list, where each list element is one variance component and if multitrait model is pursued each element of the list is a matrix of variance covariance components among traits.

iters

Maximum number of iterations allowed. Default value is 15.

tolpar

Convergence criteria.

tolparinv

tolerance parameter for matrix inverse

draw

a TRUE/FALSE value indicating if a plot showing the progress of the variance components should be drawn or not. Default is FALSE

silent

a TRUE/FALSE value indicating if the function should draw the progress bar as iterations progress. The default is FALSE, which means is not silent and will display the progress bar.

constraint

a TRUE/FALSE value indicating if the function should apply the boundary constraint indicating that variance components that are zero should be restrained.

EIGEND

a TRUE/FALSE value indicating if an eigen decomposition for the additive relationship matrix should be performed or not. This is based on Lee (2015). The limitations of this method are: 1) can only be applied to one relationship matrix 2) The system needs to be squared and no missing data is allowed (then missing data is imputed with the median). The default is FALSE to avoid the user get into trouble but experimented users can take advantage from this feature to fit big models, i.e. 5000 individuals in 555 seconds = 9 minutes in a MacBook 4GB RAM. If you use this feature keep in mind that residuals and fitted values cannot be used given the roation of the data.

forced

a list of values for variance-covariance components to be used if the user wants to force those values. The format is the same than the init argument.

complete

a TRUE/FALSE statement to indicate if the function should impute the cases where at least for one trait there's an observation. Default is TRUE in order to use all data points where at least one observation exist across traits.

restrained

a numeric argument specifying which variance-covariance parameters should be restrained. Please see vctable and vctable.help functions if you want to use it although you shouldn't since is made by the mmer2 function.

na.method.X

one of the two possible values; "include" or "exclude". If "include" is selected then the function will impute the X matrices for fixed effects with the median value. If "exclude" is selected it will get rid of all rows with missing values for the X (fixed) covariates. The default is "exclude". The "include" option should be used carefully.

na.method.Y

one of the two possible values; "include" or "exclude". If "include" is selected then the function will impute the Y matrix/dataframe with the median value. If "exclude" is selected it will get rid of missing values for the estimation of variance components. The default is "exclude".

REML

a TRUE/FALSE value to indicate if REML or ML should be used for optimization. Not functional yet. Only REML available.

Value

If all parameters are correctly indicated the program will return a list with the following information:

var.comp

a list with the values of the variance-covariance components with one list element for each random effect.

V.inv

the inverse of the phenotypic variance matrix V^- = (ZGZ+R)^-1

u.hat

a list (one element for each random effect) with a data frame for trait BLUPs.

Var.u.hat

a list (one element for each random effect) with the variance-covariance matrix for trait BLUPs.

PEV.u.hat

a list (one element for each random effect) with the predicted error variance matrix for trait BLUPs.

beta.hat

a data frame for trait BLUEs (fixed effects).

Var.beta.hat

a variance-covariance matrix for trait BLUEs

fish.inv

inverse of the Fisher's information.

residuals

Residual values e = Y - XB

cond.residuals

Conditional residual e.cond = Y - (XB + ZU)

LL

LogLikelihood

AIC

Akaike information criterion

BIC

Bayesian information criterion

X

incidence matrix for fixed effects

dimos

dimnensions for incidence matrix for random effects

sigma.scaled

scaled variance-covariance components

fitted.y

Fitted values y.hat=XB+Zu

fitted.u

Fitted values only across random effects u.hat=Zu.1+....+Zu.i

ZETA

Original incidence and variance covariance matrices used in the model fit.

K

variance-covariance matrix for random effects. If more than one random effect this will be the diagonal binding of individual K matrices.

method

The method for extimation of variance components specified by the user.

forced

a TRUE/FALSE statement indicating if user used the forced argument.

convergence

a TRUE/FALSE statement indicating if the model converged.

monitor

The values of log-likelihood and variance-covariance components across iterations during the REML estimation.

restrained

table of restrained parameters.

Details

The package has been equiped with several datasets to learn how to use the sommer package:

* HDdata and FDdata datasets have examples to fit half and full diallel designs.

* h2 to calculate heritability

* cornHybrid and Technow_data datasets to perform genomic prediction in hybrid single crosses

* wheatLines dataset to do genomic prediction in single crosses in species displaying only additive effects.

* CPdata dataset to fit genomic prediction models within a biparental population coming from 2 highly heterozygous parents including additive, dominance and epistatic effects. * PolyData to fit genomic prediction and GWAS analysis in polyploids.

* gryphondata data contains an example of an animal model including pedigree information.

* BTdata dataset contains an animal (birds) model.

If you need to use pedigree you need to convert your pedigree into a relationship matrix (use the `getA` function from the pedigreemm package). Other functions such as summary, fitted, randef (notice here is randef not ranef), anova, residuals, coef and plot applicable to typical linear models can also be applied to models fitted using the mmer-type of functions. Additional functions for genetic analysis have been included such as False Discovery Rate calculation (fdr), plot of genetic maps (map.plot), creation of manhattan plots (manhattan), phasing F1 or CP populations (phase.F1). We have also included the famous transp function to get transparent colors. Useful functions for analyzing such designs are included such as the blocker and fill.design.

MODELS ENABLED

For details about the models enabled please check the help page ot the package sommer.

References

Covarrubias-Pazaran G. Genome assisted prediction of quantitative traits using the R package sommer. PLoS ONE 2016, 11(6): doi:10.1371/journal.pone.0156744

Bernardo Rex. 2010. Breeding for quantitative traits in plants. Second edition. Stemma Press. 390 pp.

Gilmour et al. 1995. Average Information REML: An efficient algorithm for variance parameter estimation in linear mixed models. Biometrics 51(4):1440-1450.

Kang et al. 2008. Efficient control of population structure in model organism association mapping. Genetics 178:1709-1723.

Lee et al. 2015. MTG2: An efficient algorithm for multivariate linear mixed model analysis based on genomic information. Cold Spring Harbor. doi: http://dx.doi.org/10.1101/027201.

Maier et al. 2015. Joint analysis of psychiatric disorders increases accuracy of risk prediction for schizophrenia, bipolar disorder, and major depressive disorder. Am J Hum Genet; 96(2):283-294.

Searle. 1993. Applying the EM algorithm to calculating ML and REML estimates of variance components. Paper invited for the 1993 American Statistical Association Meeting, San Francisco.

Yu et al. 2006. A unified mixed-model method for association mapping that accounts for multiple levels of relatedness. Genetics 38:203-208.

Tunnicliffe W. 1989. On the use of marginal likelihood in time series model estimation. JRSS 51(1):15-27.

Zhang et al. 2010. Mixed linear model approach adapted for genome-wide association studies. Nat. Genet. 42:355-360.

Examples

# 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 for 3 variance component
#### with one variance covariance structure
####=========================================####
####=========================================####
data(CPdata)
head(CPpheno)
CPgeno[1:4,1:4]
#### 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) + Rowf + Colf, 
              rcov=~units,
              G=list(id=A), 
              data=CPpheno)
summary(mix1)
#### calculate heritability
pin(mix1, h1 ~ V1/(V1+V4) )
# assumes diag(trait) structure for univariate models

# #### for multivariate models
# mix2 <- mmer2(cbind(Yield,color)~1,
#               random = ~ us(trait):g(id) + 
#                             diag(trait):Rowf + 
#                                 diag(trait):Colf, 
#               rcov = ~ us(trait):units,
#               G=list(id=A), 
#               data=CPpheno)
# summary(mix2)

########################################################
########################################################
########################################################
########################################################
########################################################
########################################################

####=========================================####
####=========================================####
#### EXAMPLE 2
#### breeding values with 3 variance components
#### for hybrid prediction
####=========================================####
####=========================================####

# data(cornHybrid)
# hybrid2 <- cornHybrid$hybrid # extract cross data
# A <- cornHybrid$K
# K1 <- A[levels(hybrid2$GCA1), levels(hybrid2$GCA1)]; dim(K1)
# K2 <- A[levels(hybrid2$GCA2), levels(hybrid2$GCA2)]; dim(K2)
# S <- kronecker(K1, K2, make.dimnames=TRUE) ; dim(S)   
# 
# head(hybrid2)
# ans <- mmer2(Yield ~ Location, 
#              random = ~ g(GCA1) + g(GCA2) + g(SCA), 
#              rcov = ~ units,
#              G=list(GCA1=K1, GCA2=K2, SCA=S),
#              data=hybrid2)
# ans$var.comp
# summary(ans)
# 
# #### you can specify heterogeneus variances
# ans <- mmer2(Yield ~ 1, 
#              random = ~ GCA2 + at(Location):GCA2, 
#              rcov = ~ at(Location):units,
#              data=hybrid2)
# summary(ans)
# 
# ##### example of multivariate model
# ans <- mmer2(cbind(Yield,PlantHeight) ~ 1, 
#              random = ~ us(trait):GCA2 + us(trait):at(Location):GCA2, 
#              rcov = ~ diag(trait):at(Location):units,
#              data=hybrid2)
# summary(ans)

########################################################
########################################################
########################################################
########################################################
########################################################
########################################################

####=========================================####
#### EXAMPLE 3
#### different models with sommer
####=========================================####

# data(example)
# head(example)
# 
# #### Univariate homogeneous variance models ####
# 
# ## Compound simmetry (CS) model
# ans1 <- mmer2(Yield~Env, 
#               random= ~ Name + Env:Name,
#               rcov= ~ units,
#               data=example)
# summary(ans1)
# 
# #### Univariate heterogeneous variance models ####
# 
# ## Compound simmetry (CS) + Diagonal (DIAG) model
# ans2 <- mmer2(Yield~Env, 
#               random= ~Name + at(Env):Name,
#               rcov= ~ at(Env):units,
#               data=example)
# summary(ans2)
# 
# ##### Multivariate homogeneous variance models ####
# 
# ## Multivariate Compound simmetry (CS) model
# ans3 <- mmer2(cbind(Yield, Weight) ~ Env, 
#               random= ~ us(trait):Name + us(trait):Env:Name,
#               rcov= ~ us(trait):units,
#               data=example)
# summary(ans3)
# 
# ##### Multivariate heterogeneous variance models ####
# 
# ## Multivariate Compound simmetry (CS) + Diagonal (DIAG) model
# ans4 <- mmer2(cbind(Yield, Weight) ~ Env, 
#               random= ~ us(trait):Name + us(trait):at(Env):Name,
#               rcov= ~ us(trait):at(Env):units,
#               data=example)
# summary(ans4)

########################################################
########################################################
########################################################
########################################################
########################################################
########################################################

####=========================================####
#### EXAMPLE 4
#### 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,
#              rcov = ~ units,
#              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,
#               rcov = ~ at(Location):units,
#               data=hybrid2)
# summary(out2)



# }