gaston (version 1.4.9)

lmm.diago.likelihood: 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

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

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

Details

Compute the Restricted Likelihood under the linear mixed model $$ Y = (X|PC)\beta + \omega + \varepsilon $$ 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.

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")