sommer (version 1.3)

EM: Expectation Maximization Algorithm

Description

This function is used internally in the function mmer when MORE than 1 variance component needs to be estimated through the use of the expectation maximization (EM) algorithm.

Usage

EM(y, X=NULL, ETA=NULL, R=NULL, init=NULL, 
   iters=50, REML=TRUE, draw=TRUE, silent=FALSE)

Arguments

y

a numeric vector for the response variable

X

an incidence matrix for fixed effects.

ETA

an incidence matrix for random effects. This can be for one or more random effects. This NEEDS TO BE PROVIDED AS A LIST STRUCTURE. For example Z=list(list(Z=Z1, K=K1), list(Z=Z2, K=K2), list(Z=Z3, K=K3)) makes a 2 level list for 3 random effects. The gene

R

a matrix for variance-covariance structures for the residuals, i.e. for longitudinal data. if not passed is assumed an identity matrix.

init

Initial values for the variance components if defaults want to be changed.

iters

a scalar value indicating how many iterations have to be performed if the EM is performed. There is no rule of tumb for the number of iterations. The default value is 100 iterations or EM steps.

REML

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

draw

a TRUE/FALSE value indicating if a plot of updated values for the variance components should be drawn or not. The default is TRUE. COMPUTATION TIME IS SMALLER IF YOU DON'T PLOT SETTING draw=FALSE

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]

Details

This algorithm is based on Searle (1993) and Bernanrdo (2010), it is an iterative procedure NOT based on likelihood. 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 with their respective variance-covariance matrix in a 2 level list structure. For example imagine a mixed model with the following design:

.

fixed = only intercept....................................b ~ N[b.hat, 0]

random = GCA1 + GCA2 + SCA.................u ~ N[0, G]

.

where G is:

.

|K*sigma(gca1).....................0..........................0.........| |.............0.............S*sigma(gca2).....................0.........| = G

|.............0....................0......................W*sigma(sca)..|

.

The function is based on useing initial values for variance components, i.e.:

.

var(e) <- 100

var(u1) <- 100 with incidence matrix Z1

var(u2) <- 100 with incidence matrix Z2

var(u3) <- 100 with incidence matrix Z3

.

and estimates the lambda(vx) values in the mixed model equations (MME) developed by Henderson (1975), i.e. consider the 3 random effects stated above, the MME are:

.

|...............X'*R*X...............X'*R*Z1.............X'*R*Z2...................X'*R*Z3 ..............| |...X'Ry...|

|.............Z1'*R*X.........Z1'*R*Z1+K1*v1.....Z1'*R*Z2..................Z1'*R*Z3.............| |...Z1'Ry...|

|.............Z2'*R*X.............Z2'*R*Z1.............Z2'*R*Z2+K2*v2......Z2'*R*Z3.............| |...Z2'Ry...|

|.............Z3'*R*X.............Z3'*R*Z1.............Z3'*R*Z2.............Z3'*R*Z3+K3*v3......| |...Z3'Ry...|

.

..............................................................C.inv...................................................................RHS

.

where "*"" is a matrix product, R is the inverse of the var-cov matrix for the errors, Z1, Z2, Z3 are incidence matrices for random effects, X is the incidence matrix for fixed effects, K1,K2, K3 are the var-cov matrices for random effects and v1,v2,v3 are the estimates of variance components.

.

The algorithm can be summarized in the next steps:

.

1) provide initial values for the variance components

2) estimate the coefficient matrix from MME known as "C"

3) solve the mixed equations as theta = RHS * C.inv

4) obtain new estimates of fixed (b's) and random effects (u's) called theta

5) update values for variance components according to formulas

6) steps are repeated for a number of iterations specified by the user, ideally is enough when no more variations in the estimates is found, in several problems that could take thousands of iterations, whereas in other 10 iterations could be enough.

References

Bernardo Rex. 2010. Breeding for quantitative traits in plants. Second edition. Stemma Press. 390 pp.

Searle. 1993. Applying the EM algorithm to calculating ML and REML estimates of variance components. Paper invited for the 1993 American Statistical Association Meeting, San Francisco.

Examples

####=========================================####
#### For CRAN time limitations most lines in the 
#### examples are silenced with one '#' mark, 
#### remove them and run the examples
####=========================================####

###################################################
###################################################
# IMPORT DATA FOR  ESTIMATING 3 VARIANCE COMPONENTS
###################################################
###################################################
## Import phenotypic data on inbred performance
## Full data
data(cornHybrid)
hybrid2 <- cornHybrid$hybrid # extract cross data
A <- cornHybrid$K # extract the var-cov K
############################################
############################################
## breeding values with 3 variance components
############################################
############################################
y <- hybrid2$Yield
X1 <- model.matrix(~ Location, data = hybrid2);dim(X1)
Z1 <- model.matrix(~ GCA1 -1, data = hybrid2);dim(Z1)
Z2 <- model.matrix(~ GCA2 -1, data = hybrid2);dim(Z2)
Z3 <- model.matrix(~ SCA -1, data = hybrid2);dim(Z3)

K1 <- A[levels(hybrid2$GCA1), levels(hybrid2$GCA1)]; dim(K1)     
## Realized IBS relationships for set of parents 1
K2 <- A[levels(hybrid2$GCA2), levels(hybrid2$GCA2)]; dim(K2)     
## Realized IBS relationships for set of parents 2
S <- kronecker(K1, K2) ; dim(S)   
## Realized IBS relationships for cross (as the Kronecker product of K1 and K2)
rownames(S) <- colnames(S) <- levels(hybrid2$SCA)

ETA <- list(list(Z=Z1, K=K1), list(Z=Z2, K=K2), list(Z=Z3, K=S))
#ans <- EM(y=y, ETA=ETA, iters=3)