MMEst: 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)

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

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.

Author

F. Laporte and T. Mary-Huard

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.

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

Run this code

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]])

Run the code above in your browser using DataLab