mixed.solve | R Documentation |
Calculates maximum-likelihood (ML/REML) solutions for mixed models of the form
y = X \beta + Z u + \varepsilon
where \beta
is a vector of fixed effects and u
is a vector of random effects with
Var[u] = K \sigma^2_u
. The residual variance is Var[\varepsilon] = I \sigma^2_e
. This class
of mixed models, in which there is a single variance component other than the residual error,
has a close relationship with ridge regression (ridge parameter \lambda = \sigma_e^2 / \sigma^2_u
).
mixed.solve(y, Z=NULL, K=NULL, X=NULL, method="REML",
bounds=c(1e-09, 1e+09), SE=FALSE, return.Hinv=FALSE)
y |
Vector ( |
Z |
Design matrix ( |
K |
Covariance matrix ( |
X |
Design matrix ( |
method |
Specifies whether the full ("ML") or restricted ("REML") maximum-likelihood method is used. |
bounds |
Array with two elements specifying the lower and upper bound for the ridge parameter. |
SE |
If TRUE, standard errors are calculated. |
return.Hinv |
If TRUE, the function returns the inverse of |
This function can be used to predict marker effects or breeding values (see examples). The numerical method
is based on the spectral decomposition of Z K Z'
and S Z K Z' S
, where S = I - X (X' X)^{-1} X'
is
the projection operator for the nullspace of X
(Kang et al., 2008). This algorithm generates the inverse phenotypic covariance matrix V^{-1}
, which can then be used to calculate the BLUE and BLUP solutions for the fixed and random effects, respectively, using standard formulas (Searle et al. 1992):
BLUE(\beta) = \beta^* = (X'V^{-1}X)^{-1}X'V^{-1}y
BLUP(u) = u^* = \sigma^2_u KZ'V^{-1}(y-X\beta^*)
The standard errors are calculated as the square root of the diagonal elements of the following matrices (Searle et al. 1992):
Var[\beta^*] = (X'V^{-1}X)^{-1}
Var[u^*-u] = K \sigma^2_u - \sigma^4_u KZ'V^{-1}ZK + \sigma^4_u KZ'V^{-1}XVar[\beta^*]X'V^{-1}ZK
For marker effects where K = I, the function will run faster if K is not passed than if the user passes the identity matrix.
If SE=FALSE, the function returns a list containing
estimator for \sigma^2_u
estimator for \sigma^2_e
BLUE(\beta
)
BLUP(u
)
maximized log-likelihood (full or restricted, depending on method)
If SE=TRUE, the list also contains
standard error for \beta
standard error for u^*-u
If return.Hinv=TRUE, the list also contains
the inverse of H
Kang et al. 2008. Efficient control of population structure in model organism association mapping. Genetics 178:1709-1723.
Endelman, J.B. 2011. Ridge regression and other kernels for genomic selection with R package rrBLUP. Plant Genome 4:250-255.
Searle, S.R., G. Casella and C.E. McCulloch. 1992. Variance Components. John Wiley, Hoboken.
#random population of 200 lines with 1000 markers
M <- matrix(rep(0,200*1000),200,1000)
for (i in 1:200) {
M[i,] <- ifelse(runif(1000)<0.5,-1,1)
}
#random phenotypes
u <- rnorm(1000)
g <- as.vector(crossprod(t(M),u))
h2 <- 0.5 #heritability
y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g)))
#predict marker effects
ans <- mixed.solve(y,Z=M) #By default K = I
accuracy <- cor(u,ans$u)
#predict breeding values
ans <- mixed.solve(y,K=A.mat(M))
accuracy <- cor(g,ans$u)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.