In charlotte-ngs/lbgfs2020: Livestock Breeding and Genomics (Fall 2020)

knitr::opts_chunk$set(echo = FALSE, results = "asis")
knitr::knit_hooks$set(hook_convert_odg = rmdhelp::hook_convert_odg)

General Principle

• All methods to predict breeding values follow the same principle

• Correct information sources for some population mean

• Multiply corrected information source by an appropriate factor

• Regression Method $$\hat{u} = b (y - \mu)$$

• Selection Index

  • will be presented later
  • corresponds to multiple regression approach $$ \hat{u} = I = b^T y^*$$

where $b = P^{-1}Gw$ and $y^*$ corrected information sources.

Problem with Correction

• Population mean is ideal as correction

$$ y = \mu + u + e \qquad \rightarrow \qquad \bar{y} = \bar{\mu} + \bar{u} + \bar{e} = \mu $$

• Because performances are observed in different

  • environments and
  • time points

• Formation of comparison groups where animals are exposed to the same environments

• The more groups, the better the correction of environmental effects
• The more groups, the smaller the single groups

Bias

• With small comparison groups, it is more likely that mean breeding value of animals in a single group is not $0$
• Average performance of all animals in a comparison group

$$\bar{y}{CG} = \mu + \bar{u}{CG} + \bar{e}{CG}$$ * If $\bar{u}{CG}$ is not $0$, the predicted breeding value $\hat{u}_i$ of animal $i$ is

\begin{align} \hat{u}i = I &= b(y_i - (\mu + \bar{u}{CG})) \notag \ &= b(y_i - \mu) - b\bar{u}{CG} \notag \ &= \hat{u}_i - b\bar{u}{CG} \notag \end{align}

where $b\bar{u}_{CG}$ is called bias.

Solution - BLUP

• Solution to correction problem in selection index: BLUP
• Estimates environmental effects at the same time as breeding values are predicted
• Linear mixed effects model
• Meaning of BLUP
  • B stands for best $\rightarrow$ correlation between true ($u$) and its prediction ($\hat{u}$) is maximal or the prediction error variance ($var(u - \hat{u})$) is minimal.
  • L stands for linear $\rightarrow$ predicted breeding values are linear functions of the observations ($y$)
  • U stands for unbiased $\rightarrow$ expected values of the predicted breeding values are equal to the true breeding values
  • P stands for prediction

Example

### # fix the numbers parents and offspring
n_nr_sire <- 3
n_nr_dam <- 8
n_nr_parents <- n_nr_sire + n_nr_dam
n_nr_offspring <- 16
n_nr_animals <- n_nr_parents + n_nr_offspring
### # assign parents to offspring and herds to records
vec_sire_id <- c(rep(1,8), rep(2,6), rep(3,2))
vec_dam_id <- rep(4:11,2)
vec_herd_codes <- c(rep(1,4), rep(2,4), rep(1,4), rep(2,4))
### # vector of observations
vec_weaning_weight <-  c(2.61,2.31,2.44,2.41,2.51,2.55,2.14,2.61,2.34,1.99,3.1,2.81,2.14,2.41,2.54,3.16)

### # create a tibble from the data
tbl_beef_data <- dplyr::data_frame( Animal = (n_nr_parents + 1):n_nr_animals,
                                    Sire   = vec_sire_id,
                                    Dam    = vec_dam_id[order(vec_dam_id)],
                                    Herd   = vec_herd_codes,
                                    `Weaning Weight` = vec_weaning_weight )
### # count number of observations
n_nr_observation <- nrow(tbl_beef_data)

### # parameters
h2 <- .25
n_var_p <- round(var(tbl_beef_data$`Weaning Weight`), digits = 4)
n_var_g <- round(h2 * n_var_p, digits = 4)
n_pop_mean <- round(mean(tbl_beef_data$`Weaning Weight`), digits = 2)

### # show the data frame
knitr::kable( tbl_beef_data, 
              booktabs = TRUE )

Linear Models

• Simple linear model

$$y_{ij} = \mu + herd_j + e_{ij}$$

• Result: Estimate of effect of herd j
• What about breeding value $u_i$ for animal $i$?
  • Problem: breeding values have a variance $\sigma_u^2$
  • Cannot be specified in simple linear model

$\rightarrow$ Linear Mixed Effects Model (LME)

$$y_{ijk} = \mu + \beta_j + u_i + e_{ijk}$$

Matrix-Vector Notation

• LME for all animals of a population

$\rightarrow$ use matrix-vector notation

$$y = X\beta + Zu + e$$

\begin{tabular}{llp{8cm}} where & & \ & $y$ & vector of length $n$ of all observations \ & $\beta$ & vector of length $p$ of all fixed effects \ & $X$ & $n \times p$ design matrix linking the fixed effects to the observations \ & $u$ & vector of length $n_u$ of random effects \ & $Z$ & $n \times n_u$ design matrix linking random effect to the observations \ & $e$ & vector of length $n$ of random residual effects.
\end{tabular}

Expected Values and Variances

• Expected values

$$E(u) = 0 \text{ and } E(e) = 0 \rightarrow E(y) = X\beta$$

• Variances

$$var(u) = G \text{ and } var(e) = R$$

with $cov(u, e^T) = 0$, $$var(y) = Z * var(u) * Z^T + var(e) = ZGZ^T + R = V$$

The Solution

$$\hat{u} = GZ^TV^{-1}(y - X\hat{\beta})$$

$$\hat{\beta} = (X^T V^{-1} X)^- X^T V^{-1} y$$

Mixed Model Equations

$$\left[ \begin{array}{lr} X^T R^{-1} X & X^T R^{-1} Z \ Z^T R^{-1} X & Z^T R^{-1} Z + G^{-1} \end{array} \right] \left[ \begin{array}{c} \hat{\beta} \ \hat{u} \end{array} \right] = \left[ \begin{array}{c} X^T R^{-1} y \ Z^T R^{-1} y \end{array} \right]$$

Sire Model

• Breeding value of sire as random effect:

$$y = X\beta + Zs + e$$

Example

mat_x_sire <- matrix(data = c(1, 0,
                              1, 0,
                              1, 0,
                              1, 0,
                              0, 1,
                              0, 1,
                              0, 1,
                              0, 1,
                              1, 0,
                              1, 0,
                              1, 0,
                              1, 0,
                              0, 1,
                              0, 1,
                              0, 1,
                              0, 1), ncol = 2, byrow = TRUE)
vec_betahat_sire <- c("\\beta_1", "\\beta_2")
mat_z_sire <- matrix(data = c(1, 0, 0,
                              1, 0, 0,
                              1, 0, 0,
                              1, 0, 0,
                              1, 0, 0,
                              1, 0, 0,
                              1, 0, 0,
                              1, 0, 0,
                              0, 1, 0,
                              0, 1, 0,
                              0, 1, 0,
                              0, 1, 0,
                              0, 1, 0,
                              0, 1, 0,
                              0, 0, 1,
                              0, 0, 1), ncol = 3, byrow = TRUE)
vec_sirehat_sire <- c("s_1", "s_2", "s_3")
vec_res_sire <- c("e_1", "e_2", "e_3", "e_4", "e_5", "e_6", "e_7", "e_8", "e_9", "e_{10}", "e_{11}", "e_{12}", "e_{13}", "e_{14}", "e_{15}", "e_{16}")
cat("$$ \n")
# cat(paste(rmdhelp::bcolumn_vector(pvec_avector = vec_weaning_weight), sep = "\n"),"\n")
cat(paste(rmdhelp::bcolumn_vector(pvec = vec_weaning_weight), sep = "\n"), "\n")
cat("= \n")
cat(paste(rmdhelp::bmatrix(pmat = mat_x_sire), sep = "\n"), "\n")
cat(paste(rmdhelp::bcolumn_vector(pvec = vec_betahat_sire), sep = "\n"), "\n")
cat("+ \n")
cat(paste(rmdhelp::bmatrix(pmat = mat_z_sire), sep = "\n"), "\n")
cat(paste(rmdhelp::bcolumn_vector(pvec = vec_sirehat_sire), sep = "\n"), "\n")
cat("+ \n")
cat(paste(rmdhelp::bcolumn_vector(pvec = vec_res_sire), sep = "\n"), "\n")
cat("$$ \n")

Animal Model

• Breeding value for all animals as random effects

$$y = X\beta + Zu + e$$

charlotte-ngs/lbgfs2020 documentation built on Dec. 20, 2020, 5:39 p.m.