lmm.aireml | R Documentation |
Estimate the parameters of a linear mixed model, using Average Information Restricted Maximum Likelihood (AIREML) algorithm.
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)
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 |
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 |
constraint |
If |
max_iter |
Maximum number of iterations |
eps |
The algorithm stops when the gradient norm is lower than this parameter |
verbose |
If |
contrast |
If |
get.P |
If |
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.
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 |
Hervé Perdry and Claire Dandine-Roulland
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
# 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.