This function is used internally in the function `mmer`

when MORE than 1 variance component needs to be estimated through the use of the Newton-Raphson (MNR) algorithm for multivariate models (multiple responses).

1 2 3 |

`Y` |
a matrix or data frame of response variables |

`X` |
an incidence matrix for fixed effects. |

`ZETA` |
incidence matrices and var-cov matrices for random effects. This works for ONE OR MORE random effects. ' ETA <- list( A=list(Z=Z1, K=K1), B=list(Z=Z2, K=K2), C=list(Z=Z3, K=K3) ) ' makes a 2 level list for 3 the random effects A, B and C, stored in a variable we call ETA. The general idea is that each random effect is a list, i.e. A=list(Z=Z1, K=K1) where Z is the incidence matrix and K the var-cov matrix for the random effect, ' PLEASE remember to ' To provide extra detail, I'll rephrase it; when moving to situations of more than one random effect, you need to build a list for each random effect, and at the end everything gets joined in a list as well (BGLR type of format). Is called a 2-level list, i.e. A=list(Z=Z1, K=K1) and B=list(Z=Z2, K=K2) refers to 2 random effects and they should be put together in a list: ' ETA <- list( A=list(Z=Z1, K=K1), B=list(Z=Z1, K=K1) ) ' Now you can fit your model as: ' mod1 <- mmer(y=y, Z=ETA) ' You can see the examples at the bottom to have a clearer idea how to fit your models. |

`R` |
list of R matrices to correct for residual. Internally the program will do the kronecker product of such matrices to create R. |

`init` |
vector of initial values for the variance components. By default this is NULL and variance components are estimated by the method selected, but in case the user want to provide initial values this argument is functional. |

`maxcyc` |
Maximum number of cycles allowed. Default value is 50. A warning is output to the screen if this is reached before convergence. |

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

`tol` |
Convergence criteria. If the change in residual log
likelihood for one cycle is less than |

`tolparinv` |
tolerance parameter for matrix inverse |

`silent` |
a TRUE/FALSE value indicating if the function should draw the progress bar while working or should not be displayed. The default is FALSE, which means is not silent and will display the progress bar. |

`constraint` |
a TRUE/FALSE value indicating if the function should apply the boundary constraint indicating that variance components that are zero should be removed from the analysis and variance components recalculated. |

`EIGEND` |
a TRUE/FALSE value indicating if an eigen decomposition for the additive relationship matrix should be performed or not. This is based on Lee (2015). The limitations of this method are: 1) can only be applied to one relationship matrix 2) The system needs to be squared and no missing data is allowed (then missing data is imputed with the median). The default is FALSE to avoid the user get into trouble but experimented users can take advantage from this feature to fit big models, i.e. 5000 individuals in 555 seconds = 9 minutes in a MacBook 4GB RAM. |

`forced` |
a list of values for variance-covariance components to be used if the user wants to force those values. |

`IMP` |
a TRUE/FALSE statement if the function should impute the Y matrix/dataframe with the median value or should get rid of missing values for the estimation of variance components. |

This algorithm is based on Tunnicliffe (1989), it is based on REML. 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 likelihood function optimized in this algorithm is:

.

logL = -0.5 * (log( | V | ) + log( | X'VX | ) + y'Py

.

where: | | refers to the derminant of a matrix

.

If all parameters are correctly indicated the program will return a list with the following information:

- $Vu
a scalar value for the variance component estimated

- $Ve
a scalar value for the error variance estimated

- $V.inv
a matrix with the inverse of the phenotypic variance V = ZGZ+R, V^-1

- $u.hat
a vector with BLUPs for random effects

- $Var.u.hat
a vector with variances for BLUPs

- $PEV.u.hat
a vector with predicted error variance for BLUPs

- $beta.hat
a vector for BLUEs of fixed effects

- $Var.beta.hat
a vector with variances for BLUEs

- $X
incidence matrix for fixed effects, if not passed is assumed to only include the intercept

- $Z
incidence matrix for random effects, if not passed is assumed to be a diagonal matrix

- $K
the var-cov matrix for the random effect fitted in Z

- $ll
the log-likelihood value for obtained when optimizing the likelihood function when using ML or REML

Tunnicliffe W. 1989. On the use of marginal likelihood in time series model estimation. JRSS 51(1):15-27.

Covarrubias-Pazaran G (2016) Genome assisted prediction of quantitative traits using the R package sommer. PLoS ONE 11(6): doi:10.1371/journal.pone.0156744

The core functions of the package `mmer`

and `mmer2`

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 | ```
####=========================================####
#### For CRAN time limitations most lines in the
#### examples are silenced with one '#' mark,
#### remove them and run the examples
####=========================================####
data(CPdata)
CPpheno <- CPdata$pheno
CPgeno <- CPdata$geno
### look at the data
head(CPpheno)
CPgeno[1:5,1:5]
## fit a model including additive and dominance effects
Y <- CPpheno
Za <- diag(dim(Y)[1])
A <- A.mat(CPgeno) # additive relationship matrix
####================####
#### ADDITIVE MODEL ####
####================####
ETA.A <- list(add=list(Z=Za,K=A))
#ans.A <- MNR(Y=Y, ZETA=ETA.A)
#ans.A$var.comp
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

Please suggest features or report bugs with the GitHub issue tracker.

All documentation is copyright its authors; we didn't write any of that.