lmm.diago.likelihood | R Documentation |
Compute the Restricted or the Full Likelihood of a linear mixed model, using the "diagonalization trick".
lmm.diago.likelihood(tau, s2, h2, Y, X, eigenK, p = 0) lmm.diago.profile.likelihood(tau, s2, h2, Y, X, eigenK, p = 0)
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 |
Theses function respectively compute the Restricted and the Profile 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, 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/(tau + sigma^2) = h^2,
and output these values as well as the likelihood value at this point.
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 |
Hervé Perdry and Claire Dandine-Roulland
lmm.restricted.likelihood
, lmm.profile.restricted.likelihood
, lmm.diago
, lmm.aireml
# 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")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.