# lmm.aireml: Linear mixed model fitting with AIREML In gaston: Genetic Data Handling (QC, GRM, LD, PCA) & Linear Mixed Models

## Description

Estimate the parameters of a linear mixed model, using Average Information Restricted Maximum Likelihood (AIREML) algorithm.

## Usage

 1 2 3 4 5 lmm.aireml(Y, X = matrix(1, nrow = length(Y)), K, EMsteps = 0L, EMsteps_fail = 1L, EM_alpha = 1, min_tau, min_s2 = 1e-06, theta, constraint = TRUE, max_iter = 50L, eps = 1e-05, verbose = getOption("gaston.verbose", TRUE), contrast = FALSE, get.P = FALSE) 

## Arguments

 Y Phenotype vector X Covariable matrix. By default, a column of ones to include an intercept in the model K A positive definite matrix or a list of such matrices EMsteps Number of EM steps ran prior the AIREML EMsteps_fail Number of EM steps performed when the AIREML algorithm fail to improve the likelihood value EM_alpha Tweaking parameter for the EM (see Details) min_tau Minimal value for model parameter tau (if missing, will be set to 1e-6) min_s2 Minimal value for model parameter sigma^2 theta (Optional) Optimization starting point theta = c(sigma^2, tau) constraint If TRUE, the model parameters respect the contraints given by min_tau and min_s2 max_iter Maximum number of iterations eps The algorithm stops when the gradient norm is lower than this parameter verbose If TRUE, display information on the algorithm progress contrast If TRUE, use a contrast matrix to compute the Restricted Likelihood (usually slower) get.P If TRUE, the function sends back the last matrix P computed in the optimization process

## Details

Estimate the parameters of the following linear mixed model, using AIREML algorithm:

Y =X beta + omega_1 + ... + omega_k + epsilon

with omega_i ~ N(0, tau_i K_i) for i \in 1, …,k and epsilon ~ N(0, sigma^2 I_n).

The variance matrices K_1, ..., K_k, are specified through the parameter K.

If EMsteps is positive, the function will use this number of EM steps to compute a better starting point for the AIREML algorithm. Setting EMsteps to a value higher than max_iter leads to an EM optimization. It can happen that after an AIREML step, the likelihood did not increase: if this happens, the functions falls back to EMsteps_fail EM steps. The parameter EM_alpha can be set to a value higher than 1 to attempt to accelerate EM convergence; this could also result in uncontrolled behaviour and should be used with care.

After convergence, the function also compute Best Linear Unbiased Predictors (BLUPs) for beta and omega, and an estimation of the participation of the fixed effects to the variance of Y.

## Value

A named list with members:

 sigma2 Estimate of the model parameter sigma^2 tau Estimate(s) of the model parameter(s) tau_1, ..., tau_k logL Value of log-likelihood logL0 Value of log-likelihood under the null model (without random effect) niter Number of iterations done norm_grad Last computed gradient's norm P Last computed value of matrix P (see reference) Py Last computed value of vector Py (see reference) BLUP_omega BLUPs of random effects BLUP_beta BLUPs of fixed effects beta varbeta Variance matrix for beta estimates varXbeta Participation of fixed effects to variance of Y

If get.P = TRUE, there is an additional member:

 P The last matrix P computed in the AIREML step

## Author(s)

Herv<c3><a9> Perdry and Claire Dandine-Roulland

## References

Gilmour, A. R., Thompson, R., & Cullis, B. R. (1995), Average information REML: an efficient algorithm for variance parameter estimation in linear mixed models, Biometrics, 1440-1450

lmm.diago, logistic.mm.aireml, lmm.simu
  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 # Temporarily set nb threads to 2 to comply with CRAN rules n.threads <- getNumThreads() setThreadOptions(2) # Load data data(AGT) x <- as.bed.matrix(AGT.gen, AGT.fam, AGT.bim) # Compute Genetic Relationship Matrix standardize(x) <- "p" K <- GRM(x) # Simulate a quantitative genotype under the LMM set.seed(1) y <- 1 + x %*% rnorm(ncol(x), sd = 1)/sqrt(ncol(x)) + rnorm(nrow(x), sd = sqrt(2)) # Estimates estimates <- lmm.aireml(y, K = K, verbose = FALSE) str(estimates) # back to previous nb threads value setThreadOptions(n.threads)