Data Engineering and BI courses are free this week!
mmer when MORE than 1 variance component needs to be estimated through the use of the expectation maximization (EM) algorithm.EM(y, X=NULL, ETA=NULL, R=NULL, init=NULL,
iters=50, REML=TRUE, draw=TRUE, silent=FALSE, forced=NULL).
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.
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.
####=========================================####
#### 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)Run the code above in your browser using DataLab