lmm.diago.likelihood {gaston}

Likelihood of a linear mixed model

Description

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

Usage

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

Arguments

tau

Value(s) of model parameter (see Details)

s2

Value(s) of model parameter (see Details)

h2

Value(s) of heritability (see Details)

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

Details

Theses function respectively compute the Restricted and the Profile Likelihood under the linear mixed model

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.

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

The function lmm.diago.profile.likelihood computes the full likelihood, profiled for \beta. That is, the value \beta which maximizes the full likelihood for the given values of \tau and \sigma^2 is computed, and then the full likelihood is computed.

If h2 is provided, both functions compute \tau and \sigma^2 which maximizes the likelihood under the constraint {\tau \over \tau + \sigma^2 } = h^2 , and output 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

tau

Corresponding values of \tau

sigma2

Corresponding values of \sigma^2

likelihood

Corresponding likelihood values

Author(s)

Hervé Perdry and Claire Dandine-Roulland

See Also

lmm.restricted.likelihood, lmm.profile.restricted.likelihood, 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")