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 |
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 + \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.
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 |
sigma2 |
Corresponding values of |
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.