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

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 memberstauCorresponding values of $\tau$sigma2Corresponding values of $\sigma^2$likelihoodCorresponding likelihood values

Details

Compute the Restricted 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, 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 $\sigma^2$ which maximizes the likelihood under the constraint ${\tau \over \tau + \sigma^2 } = h^2$, and outputs these values as well as the likelihood value at this point.

Examples

Run this code

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