sommer (version 1.3)

EMMA: Efficient Mixed Model Association Algorithm

Description

This function is used internally in the function mmer when ONLY ONE variance component other than the error needs to be estimated through the use of the efficient mixed model association (EMMA) algorithm.

Usage

EMMA(y, X=NULL, Z=NULL, K=NULL, REML=TRUE, silent=FALSE)

Arguments

y

a numeric vector for the response variable

X

an incidence matrix for fixed effects.

Z

an incidence matrix for the random effect fitted other than the error variance component.

K

a var-cov matrix for the random effect fitted in Z.

REML

a TRUE/FALSE value indicating if restricted maximum likelihood should be used instead of ML. The default is TRUE.

silent

a TRUE/FALSE value indicating if the function should draw the progress bar or iterations performed while working or should not be displayed.

Value

• If all parameters are correctly indicated the program will return a list with the following information: [object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object]

Details

This algorithm is based on Kang et al. (2008), it is based on REML using the ridge regression parameter lambda as the ration of Var(e)/Var(u). This handles models of the form:

.

y = Xb + Zu + e

.

b ~ N[b.hat, 0] zero variance because is a fixed term

u ~ N[0, K*sigma(u)] where: K*sigma(u) = G

e ~ N[0, I*sigma(e)] where: I*sigma(e) = R

y ~ N[Xb, var(Zu+e)] where;

var(y) = var(Zu+e) = ZGZ+R = V which is the phenotypic variance

.

The function allows the user to specify the incidence matrices for the random effect "u" as Z and its variance-covariance matrix as K, only one Z and K for one random effect.

.

The likelihood function optimized in this algorithm is:

.

logL = (n - p) * log(sum(eta^2/{ lambda + delta})) + sum(log(lambda + delta))

.

where: (n-p) refers to the degrees of freedom lambda are the eigenvalues mentioned by Kang et al.(2008) delta is the REML estimator of the ridge parameter

.

The algorithm can be summarized in the next steps:

.

1) provide initial value for the ridge parameter

2) estimate S = I - X(X'X)-X'

3) obtain the phenotypic variance V = ZKZ' + delta.prov*I

4) perform an eigen decomposition of SVS

5) create "lambda"" as the eigenvalues of SVS and "U"" as the eigenvectors

6) estimate eta=U'y

7) optimize the likelihood shown above providing "eta", "lambdas" and optimize with respect to "delta" which is the ridge parameter and contains Ve/Vu

References

Kang et al. 2008. Efficient control of population structure in model organism association mapping. Genetics 178:1709-1723.

Examples

####=========================================####
#### breeding values with 1 variance component
#### using EMMA algorithm
####=========================================####

####=========================================####
#### 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)
}
####=========================================####
#### Simulate 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)))

ans <- EMMA(y=y, K=A.mat(M))