Likelihood of a linear mixed model

Description

Compute the Restricted Likelihood of a linear mixed model, using the "diagonalization trick".

Usage

 lmm.diago.likelihood(tau, s2, h2, Y, X, eigenK, p = 0) 

Arguments

Details

Compute the Restricted Likelihood under the linear mixed model

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.

If both tau and s2 (for sigma^2) are provided, the function computes the likelihood for these values of the parameters; if these parameters are vectors of length > 1, then a matrix of likelihood values is computed.

If h2 is provided, the function computes tau and s2 which maximizes the likelihood under the constraint tau/(tau + s2) = h2, and outputs these values as well as the likelihood value at this point.

Value

If tau and s2 are provided, the corresponding likelihood values.

If tau or s2 are missing, and h2 is provided, a named list with members

Author(s)

Herv<c3><83><c2><a9> Perdry and Claire Dandine-Roulland

See Also

lmm.diago, lmm.aireml

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
     
# Likelihood
TAU <- seq(0.5,1.5,length=30)
S2 <- seq(1,3,length=30)
lik1 <- lmm.diago.likelihood(tau = TAU, s2 = S2, Y = y, eigenK = eiK)

H2 <- seq(0,1,length=51)
lik2 <- lmm.diago.likelihood(h2 = H2, Y = y, eigenK = eiK)

# Plotting
par(mfrow=c(1,2))
lik.contour(TAU, S2, lik1, heat = TRUE, xlab = "tau", ylab = "sigma^2")
lines(lik2$tau, lik2$sigma2)
plot(H2, exp(lik2$likelihood), type="l", xlab="h^2", ylab = "likelihood")

Results