lmm.diago: Linear mixed model fitting with the diagonalization trick

In gaston: Genetic Data Handling (QC, GRM, LD, PCA) & Linear Mixed Models

Linear mixed model fitting with the diagonalization trick

Description

Estimate the parameters of a linear mixed model, using the "diagonalization trick".

Usage

 lmm.diago(Y, X = matrix(1, nrow=length(Y)), eigenK, p = 0, 
                  method = c("newton", "brent"), min_h2 = 0, max_h2 = 1,
                  verbose = getOption("gaston.verbose", TRUE), 
                  tol = .Machine$double.eps^0.25) 

Arguments

Y	Phenotype vector
X	Covariable matrix
eigenK	Eigen decomposition of K (a positive symmetric matrix)
p	Number of Principal Components included in the mixed model with fixed effect
method	Optimization method to use
min_h2	Minimum admissible value
max_h2	Maximum admissible value
verbose	If TRUE, display information on the function actions
tol	Accuracy of estimation

Details

Estimate the parameters of the following linear mixed model, computing the restricted likelihood as in lmm.diago.likelihood, and using either a Newton algorithm, or Brent algorithm as in optimize:

Y = (X|PC) beta + omega + epsilon

with omega ~ N(0, tau K) and epsilon ~ N(0, sigma^2 I_n).

The matrix K is given through its eigen decomposition, as produced by eigenK = eigen(K, symmetric = TRUE). The matrix (X|PC) is the concatenation of the covariable matrix X and of the first p eigenvectors of K, included in the model with fixed effects.

Value

If the parameter p is a scalar, a list with following elements :

sigma2	Estimate of the model parameter sigma^2
tau	Estimate(s) of the model parameter(s) tau_1, ..., tau_k
Py	Last computed value of vector Py (see reference)
BLUP_omega	BLUPs of random effects
BLUP_beta	BLUPs of fixed effects beta (only the components corresponding to X)
Xbeta	Estimate of (X|PC)beta
varbeta	Variance matrix for beta estimates (only the components corresponding to X)
varXbeta	Participation of fixed effects to variance of Y
p	Number of Principal Components included in the linear mixed model with fixed effect

If the paramer p is a vector of length > 1, a list of lists as described above, one for each value in p.

Author(s)

Hervé Perdry and Claire Dandine-Roulland

See Also

lmm.diago.likelihood, lmm.aireml, optimize

Examples

# Load data
data(AGT)
x <- as.bed.matrix(AGT.gen, AGT.fam, AGT.bim)

# Compute Genetic Relationship Matrix
K <- GRM(x)

# eigen decomposition of K
eiK <- eigen(K)

# simulate a phenotype
set.seed(1)
y <- 1 + lmm.simu(tau = 1, sigma2 = 2, eigenK = eiK)$y

# Estimations
R <- lmm.diago(Y = y, eigenK = eiK, p = c(0,10))
str(R)

gaston documentation built on Dec. 28, 2022, 1:30 a.m.