gaston (version 1.0)

lmm.diago: 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, 
                  tol = .Machine$double.eps^0.25)

Value

If the parameter p is a scalar, a list with following elements :sigma2Estimate of the model parameter $\sigma^2$tauEstimate(s) of the model parameter(s) $\tau_1, \dots, \tau_k$PyLast computed value of vector Py (see reference)BLUP_omegaBLUPs of random effectsBLUP_betaBLUPs of fixed effects $\beta$ (only the components corresponding to $X$)XbetaEstimate of $(X|PC)\beta$varbetaVariance matrix for $\beta$ estimates (only the components corresponding to $X$)varXbetaParticipation of fixed effects to variance of YpNumber of Principal Components included in the linear mixed model with fixed effectIf the paramer p is a vector of length > 1, a list of lists as described above, one for each value in p.

Details

Estimate the parameters of the following linear mixed model, computing the likelihood as in lmm.diago.likelihood, and using the optimization algorithm of the R base function optimize: $$Y = (X|PC)\beta + \omega + \varepsilon$$ with $\omega \sim N(0,\tau K)$ and $\varepsilon \sim 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.

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
x <- set.stats(x)
standardize(x) <- 'mu'
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)