gaston (version 1.4.9)

lmm.aireml: Linear mixed model fitting with AIREML

Description

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

Usage

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

Value

A named list with members:If get.P = TRUE, there is an additional member:

Details

Estimate the parameters of the following linear mixed model, using AIREML algorithm: $$ Y = X\beta + \omega_1 + \ldots + \omega_k + \varepsilon $$ with $omega_i ~ N(0, tau_i K_i)$ for $ i \in 1, \dots,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$.

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

See Also

lmm.diago, lmm.simu

Examples

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