MMEst {MM4LMM}

MM inference method for variance component mixed models

Description

This is the main function of the MM4LMM package. It performs inference in a variance component mixed model using a Min-Max algorithm. Inference in multiple models (e.g. for GWAS analysis) can also be performed.

Usage

MMEst(Y, Cofactor = NULL, X = NULL, formula=NULL, VarList, ZList = NULL, Method = "Reml",
	Henderson=NULL, Init = NULL, CritVar = 0.001, CritLogLik = 0.001,
	MaxIter = 100, NbCores = 1, Verbose = TRUE)

Arguments

Y

A vector of response values.

Cofactor

An incidence matrix corresponding to fixed effects common to all models to be adjusted. If NULL, a single intercept is used as cofactor.

X

An incidence matrix or a list of incidence matrices corresponding to fixed effects specific to each model. If X is a matrix, one model per column will be fitted. If X is a list, one model per element of the list will be fitted (default is NULL).

formula

A formula object specifying the fixed effect part of all models separated by + operators. To specify an interaction between Cofactor and X use the colnames of X when it is a list or use "Xeffect" when X is a matrix.

VarList

A list of covariance matrices associated with random and residual effects.

ZList

A list of incidence matrices associated with random and residual effects (default is NULL).

Method

The method used for inference. Available methods are "Reml" (Restricted Maximum Likelihood) and "ML" (Maximum Likelihood).

Henderson

If TRUE the Henderson trick is applied. If FALSE the Henderson trick is not applied. If NULL the algorithm chooses wether to use the trick or not.

Init

A vector of initial values for variance parameters (default is NULL)

CritVar

Value of the criterion for the variance components to stop iteration. (see Details)

CritLogLik

Value of the criterion for the log-likelihood to stop iteration. (see Details)

MaxIter

Maximum number of iterations per model.

NbCores

Number of cores to be used.

Verbose

A boolean describing if messages have to be printed (TRUE) or not (FALSE). Default is TRUE.

Details

If X is NULL, the following model is fitted:

Y = X_C \beta_C + \sum_{k=1}^K Z_k u_k

with X_C the matrix provided in Cofactor, \beta_C the unknown fixed effects, Z_k the incidence matrix provided for the kth component of ZList and u_k the kth vector of random effects. If ZList is unspecified, all incidence matrices are assumed to be the Identity matrix. Random effects are assumed to follow a Gaussian distribution with mean 0 and covariance matrix R_k \sigma_k^2, where R_k is the kth correlation matrix provided in VarList.

If X is not NULL, the following model is fitted for each i:

Y = X_C \beta_C + X_{[i]} \beta_{[i]} + \sum_{k=1}^K Z_k u_k

where X_{[i]} is the incidence matrix corresponding to the ith component (i.e. column if X is a matrix, element otherwise) of X, and \beta_{[i]} is the (unknow) fixed effect associated to X_{[i]}.

All models are fitted using the MM algorithm. If Henderson=TRUE, at each step the quantities required for updating the variance components are computed using the Mixed Model Equation (MME) trick. See Johnson et al. (1995) for details.

Value

The result is a list where each element corresponds to a fitted model. Each element displays the following:

Beta

Estimated values of \beta_C and \beta_{i}

Sigma2

Estimated values of \sigma_1^2,...,\sigma_K^2

VarBeta

Variance matrix of Beta

LogLik (Method)

The value of the (restricted, if Method is "Reml") log-likelihood

NbIt

The number of iterations required to reach the optimum

Method

The method used for the inference

attr

An integer vector with an entry for each element of Beta giving the term in Factors which gave rise to this element (for internal use in the function AnovaTest)

Factors

Names of each term in the formula

Author(s)

F. Laporte and T. Mary-Huard

References

Laporte, F., Charcosset, A., & Mary-Huard, T. (2022). Efficient ReML inference in variance component mixed models using a Min-Max algorithm. PLOS Computational Biology, 18(1), e1009659.

Johnson, D. L., & Thompson, R. (1995). Restricted maximum likelihood estimation of variance components for univariate animal models using sparse matrix techniques and average information. Journal of dairy science, 78(2), 449-456.

Hunter, D. R., & Lange, K. (2004). A tutorial on MM algorithms. The American Statistician, 58(1), 30-37.

Zhou, H., Hu, L., Zhou, J., & Lange, K. (2015). MM algorithms for variance components models. arXiv preprint arXiv:1509.07426.

Examples

require('MM4LMM')

#### Example 1: variance component analysis, 1 model
data(VarianceComponentExample)
DataHybrid <- VarianceComponentExample$Data
KinF <- VarianceComponentExample$KinshipF
KinD <- VarianceComponentExample$KinshipD

##Build incidence matrix for each random effect
Zf <- t(sapply(as.character(DataHybrid$CodeFlint), function(x)
  as.numeric(rownames(KinF)==x)))
Zd <- t(sapply(as.character(DataHybrid$CodeDent), function(x)
  as.numeric(rownames(KinD)==x)))

##Build the VarList and ZList objects
VL = list(Flint=KinF , Dent=KinD , Error = diag(1,nrow(DataHybrid)))
ZL <- list(Flint=Zf , Dent=Zd , Error = diag(1,nrow(DataHybrid)))

##Perform inference
#A first way to call MMEst
ResultVA <- MMEst(Y=DataHybrid$Trait  , Cofactor = matrix(DataHybrid$Trial)
                  , ZList = ZL  ,  VarList = VL)
length(ResultVA)
print(ResultVA)

#A second way to call MMEst (same result)
Formula <- as.formula('~ Trial')
ResultVA2 <- MMEst(Y=DataHybrid$Trait  , Cofactor = DataHybrid,
                   formula = Formula
                  , ZList = ZL  ,  VarList = VL)
length(ResultVA2)
print(ResultVA2)



#### Example 2: Marker Selection with interaction between Cofactor and X matrix
Formula <- as.formula('~ Trial+Xeffect+Xeffect:Trial')
ResultVA3 <- MMEst(Y=DataHybrid$Trait  , Cofactor = DataHybrid,
                  X = VarianceComponentExample$Markers,
                   formula = Formula
                  , ZList = ZL  ,  VarList = VL)
length(ResultVA3)
print(ResultVA3[[1]])


#### Example 3: QTL detection with two variance components
data(QTLDetectionExample)
Pheno <- QTLDetectionExample$Phenotype
Geno <- QTLDetectionExample$Genotype
Kinship <- QTLDetectionExample$Kinship

##Build the VarList object
VLgd <- list(Additive=Kinship , Error=diag(1,length(Pheno)))

##Perform inference
ResultGD <- MMEst(Y=Pheno , X=Geno
                  , VarList=VLgd , CritVar = 10e-5)

length(ResultGD)
print(ResultGD[[1]])