gaston (version 1.2)

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)

Value

A named list with members:sigma2Estimate of the model parameter $\sigma^2$tauEstimate(s) of the model parameter(s) $\tau_1, \dots, \tau_k$logLValue of log-likelihoodlogL0Value of log-likelihood under the null model (without random effect)niterNumber of iterations donenorm_gradLast computed gradient's normPyLast computed value of vector Py (see reference)BLUP_omegaBLUPs of random effectsBLUP_betaBLUPs of fixed effects $\beta$varbetaVariance matrix for $\beta$ estimatesvarXbetaParticipation of fixed effects to variance of Y

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 \sim N(0,\tau_i K_i)$ for $i \in 1, \dots,k$ and $\varepsilon \sim 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
x <- set.stats(x)
standardize(x) <- 'mu'
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)