Extracting the parameters and predictions from the lme4breeding models

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

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.

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)

## boundary (singular) fit: see help('isSingular')

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.

vc <- VarCorr(ans1); print(vc,comp=c("Variance"))

##  Groups    Name        Variance
##  Env:Name  (Intercept)  5.1528 
##  Name      (Intercept)  3.7184 
##  Env:Block (Intercept)  0.0000 
##  Env       (Intercept) 12.0084 
##  Residual               4.3661

ve <- attr(VarCorr(ans1), "sc")^2; ve

## [1] 4.366116

3) Extracting random effect coefficients and their standard errors

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

##            (Intercept)
## A01143-3C    2.9601577
## AC00206-2W  -0.2291471
## AC01151-5W   0.7260436
## AC03433-1W  -0.9611183
## AC03452-2W   1.4069492
## AC05153-1W  -1.4876828

head(PEV$Name)

##            (Intercept)
## A01143-3C     1.537279
## AC00206-2W    1.894723
## AC01151-5W    1.903307
## AC03433-1W    1.537279
## AC03452-2W    1.894723
## AC05153-1W    1.894723

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.

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.

pp <- predict(ans1, classify="Name")

## hyperTable argument not provided. Building a hyper table based on the classify argument. Please check the output slot 'hyperTable' to ensure that the different effects have been included and average as you expected.

head(pp$pvals)

##                    id predicted.value std.error
## A01143-3C   A01143-3C       19.488758  3.308065
## AC00206-2W AC00206-2W       11.879814  2.797327
## AC01151-5W AC01151-5W       14.158678  2.806046
## AC03433-1W AC03433-1W       10.133500  3.308065
## AC03452-2W AC03452-2W       15.783162  2.797327
## AC05153-1W AC05153-1W        8.877238  2.797327

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

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:

pp$hyperTable

##        variable       group   type include average
## 1   (Intercept) (Intercept)  fixed       1       1
## 2   (Intercept)    Env:Name random       1       1
## 96  (Intercept)        Name random       1       0
## 137 (Intercept)   Env:Block random       0       0
## 143 (Intercept)         Env random       0       0

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.