lmm.aireml {gaston} | R Documentation |
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 |
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 |
min_s2 |
Minimal value for model parameter |
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 |
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
.
Value
A named list with members:
sigma2 |
Estimate of the model parameter |
tau |
Estimate(s) of the model parameter(s) |
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 |
varbeta |
Variance matrix for |
varXbeta |
Participation of fixed effects to variance of Y |
If get.P = TRUE
, there is an additional member:
P |
The last matrix |
Author(s)
Hervé 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
See Also
lmm.diago
, logistic.mm.aireml
, 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)