---
title: "Extracting the parameters and predictions from the lme4breeding models"
author: "Giovanny Covarrubias-Pazaran"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Extracting parameters from the lme4breeding models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

lme4breeding is nice wrapper of the lme4 package that enables the use of especialized plant and animal breeding models that include relationship matrices among individuals (e.g., genomic relationship matrices) and complex covariance structures between factors (e.g., factor analytic structures). It uses all the lme4 machinery for linear and non-linear models, for different response distributions opening a world of possibilities. 

This vignette aim to help users to understand how to retrieve important parameters of relevance: 


1) Extracting fixed effect coefficients and their standard errors
2) Extracting variance-covariance components
3) Extracting random effect coefficients and their standard errors
4) Computing predictions or linear combinations

### 1) Extracting fixed effect coefficients and their standard errors

Of the basic output of a linear model is the fixed coefficients that many times are used for hypothesis testing. Here we show the use of the `fixef()` and `vcov()` functions.

```{r setup, include=FALSE} 
# knitr::opts_chunk$set(warning = FALSE, message = FALSE) 
library(lme4breeding)
```

```{r}
data(DT_example)
DT <- DT_example
A <- A_example

ans1 <- lmebreed(Yield~ (1|Name) + (1|Env) + 
                   (1|Env:Name) + (1|Env:Block),
             verbose = FALSE, data=DT)
```

### 2) Extracting variance-covariance components

To extract the variance-covariance components from a model we can use the `VarCorr` function. To extract the residual variance we need to access the attributes of a `VarCorr` object.

```{r}
vc <- VarCorr(ans1); print(vc,comp=c("Variance"))
ve <- attr(VarCorr(ans1), "sc")^2; ve
```

### 3) Extracting random effect coefficients and their standard errors

```{r}
BLUP <- ranef(ans1, condVar=TRUE)
PEV <- lapply(BLUP, function(x){attr(x, which="postVar")}) # take sqrt() for SEs
head(BLUP$Name)
head(PEV$Name)
```

Alternatively, you can extract the predicted error variance matrix (inverse of the coefficient matrix) and extract the standard error from coefficients taking the square root from the diagonal elements of the matrix.

```{r}
Ci <- getCi(ans1)
Matrix::image(Ci)
```

This matrix is also included in the attributes of the object returned by the `ranef()` function.

### 4) Computing predictions or linear combinations

One of the most searched functions in genetic evaluation is the `predict()` function which aims to compute linear combinations for fixed and random effects.

```{r}
pp <- predict(ans1, classify="Name")
head(pp$pvals)
```

The object coming out of `predict()` includes the linear combination matrix (D)

```{r}
image(pp$D)
```

and the hypertable used to create the D matrix. The user can provide their own hypertable to control which effects should be included, excluded or included and averaged:

```{r}
pp$hyperTable
```

## Literature

Giovanny Covarrubias-Pazaran (2024).  lme4breeding: enabling genetic evaluation in the age of genomic data. To be submitted.

Bates Douglas, Maechler Martin, Bolker Ben, Walker Steve. 2015. Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1-48.

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

Henderson C.R. 1975. Best Linear Unbiased Estimation and Prediction under a Selection Model. Biometrics vol. 31(2):423-447.

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.

Welham, S., Cullis, B., Gogel, B., Gilmour, A., and Thompson, R. (2004). Prediction in linear mixed models. Australian and New Zealand Journal of Statistics, 46, 325 - 347.

