In charlotte-ngs/lbgfs2020: Livestock Breeding and Genomics (Fall 2020)
knitr::opts_chunk$set(echo = FALSE, results = 'asis', fig.pos = 'H')
knitr::knit_hooks$set(hook_convert_odg = rmdhelp::hook_convert_odg)
Problem 1: Numerator Relationship Matrix
Construct the numerator relationship matrix $A$ for the following pedigree and verify the result using the function getA() from package pedigreemm.
tbl_ped_ex8p1 <- tibble::tibble(Animal = c(5, 6, 7, 8, 9, 10),
Sire = c(1, 1, 4, 4, 4, 4),
Dam = c(2, 3, 5, 5, 6, 6))
knitr::kable(tbl_ped_ex8p1,
booktabs = TRUE,
longtable = TRUE,
caption = "Pedigree For Constructing Numerator Relationship Matrix")
Solution
The numerator relationship is constructed using the follownig step-wise procedure. The following rules are used to compute the single elements.
-
Case 1: If both parents $s$ and $d$ of animal $i$ are known then
- the diagonal element $(A){ii}$ corresponds to: $(A){ii} = 1 + F_i = 1 + {1\over 2} (A)_{sd}$ and
- the offdiagonal element $(A){ji}$ is computed as: $(A){ji} = {1\over 2} ((A){js} + (A){jd})$
- because $A$ is symmetric $(A){ji} = (A){ij}$
-
Case 2: If only one parent $s$ is known and assumed unrelated to the mate
- $(A)_{ii} = 1$
- $(A){ij} = (A){ji} = {1\over 2} ((A)_{js}$
-
Case 3: If both parents are unknown
- $(A)_{ii} = 1$
- $(A){ij} = (A){ji} = 0$
Step 1
First, we extend the pedigree given in Table \@ref(tab:ex8p1ped). All animals without parents are added at the top of the pedigree. This results in the matrix shown in Table \@ref(tab:ex8p1pedext).
n_nr_ani_ex8p1pedext <- max(tbl_ped_ex8p1$Animal)
n_nr_founder_ex8p1pedext <- min(tbl_ped_ex8p1$Animal) - 1
tbl_ped_exp8p1ext <- tibble::tibble(Animal = c(1:n_nr_ani_ex8p1pedext),
Sire = c(rep(NA, n_nr_founder_ex8p1pedext), tbl_ped_ex8p1$Sire),
Dam = c(rep(NA, n_nr_founder_ex8p1pedext), tbl_ped_ex8p1$Dam))
knitr::kable(tbl_ped_exp8p1ext,
booktabs = TRUE,
longtable = TRUE,
caption = "Extended Pedigree")
Because the pedigree in Table \@ref(tab:ex8p1pedext) is already ordered such that parents are before offspring, we can directly go to the next step.
Step 2
We start with an empty numerator relationship matrix $A$. The matrix $A$ has dimensions $r n_nr_ani_ex8p1pedext \times r n_nr_ani_ex8p1pedext$
matA_empty_ex8p1 <- matrix(NA, nrow = n_nr_ani_ex8p1pedext, ncol = n_nr_ani_ex8p1pedext)
matA_empty_string_ex8p1 <- matrix('', nrow = n_nr_ani_ex8p1pedext, ncol = n_nr_ani_ex8p1pedext)
### # display
cat(paste(rmdhelp::bmatrix(pmat = matA_empty_string_ex8p1, ps_name = 'A', ps_env = '$$')))
Step 3
The single elements of $A$ are computed according to the rules listed above.
n_ex8p1_first_animal <- tbl_ped_exp8p1ext$Animal[1]
The computation is started with animal $r n_ex8p1_first_animal$. The first element is always the diagonal-element that corresponds to animal that we are currently looking at. For animal $r n_ex8p1_first_animal$ the diagonal element is $(A)_{r n_ex8p1_first_animal``r n_ex8p1_first_animal}$. Because animal r n_ex8p1_first_animal has not parents, we are in case 3 for the diagonal element. If an animal has unknown parents, it also means that the animals's inbreeding coefficient $F_i$ is $0$. Hence
$$(A)_{r n_ex8p1_first_animal``r n_ex8p1_first_animal} = 1$$
Now we have the first element of our numerator relationship matrix.
matA_empty_ex8p1[n_ex8p1_first_animal,n_ex8p1_first_animal] <- 1
matA_empty_string_ex8p1[n_ex8p1_first_animal,n_ex8p1_first_animal] <- 1
### # display
cat(paste0(rmdhelp::bmatrix(pmat = matA_empty_string_ex8p1, ps_name = 'A', ps_env = '$$')))
Step 4
The next elements that need to be computed are the off-diagnoal element on row $1$. Elements $A_{12}$, $A_{13}$ and $A_{14}$ correspond to the additive genetic relationship between animal $1$ and animals $2$, $3$ and $4$. Because animals $2$, $3$ and $4$ all have unknown parents, we are for all three elements in case 3, hence we can state
$$(A){12} = (A){13} = (A)_{14} = 0$$
matA_empty_ex8p1[n_ex8p1_first_animal,2:4] <- 0
For the remaining elements of the first row of $A$, the elements correspond to the additive genetic relationship between animal $1$ and animals $5$ to $10$. Because animals $5$ to $10$ all have known parents, we have to use case 1 in the above formulated rules.
n_cur_ani <- 5
n_cur_sire <- tbl_ped_exp8p1ext$Sire[n_cur_ani]
n_cur_dam <- tbl_ped_exp8p1ext$Dam[n_cur_ani]
matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] <- 0.5 * (matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam])
\begin{equation}
(A){r n_ex8p1_first_animal``r n_cur_ani} = {1\over 2}\left( (A){r n_ex8p1_first_animal``r n_cur_sire} + (A)_{r n_ex8p1_first_animal``r n_cur_dam} \right)
= {1\over 2}\left( r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam] \right)
= r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] \notag
\end{equation}
n_cur_ani <- 6
n_cur_sire <- tbl_ped_exp8p1ext$Sire[n_cur_ani]
n_cur_dam <- tbl_ped_exp8p1ext$Dam[n_cur_ani]
matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] <- 0.5 * (matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam])
\begin{equation}
(A){r n_ex8p1_first_animal``r n_cur_ani} = {1\over 2}\left( (A){r n_ex8p1_first_animal``r n_cur_sire} + (A)_{r n_ex8p1_first_animal``r n_cur_dam} \right)
= {1\over 2}\left( r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam] \right)
= r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] \notag
\end{equation}
n_cur_ani <- 7
n_cur_sire <- tbl_ped_exp8p1ext$Sire[n_cur_ani]
n_cur_dam <- tbl_ped_exp8p1ext$Dam[n_cur_ani]
matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] <- 0.5 * (matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam])
\begin{equation}
(A){r n_ex8p1_first_animal``r n_cur_ani} = {1\over 2}\left( (A){r n_ex8p1_first_animal``r n_cur_sire} + (A)_{r n_ex8p1_first_animal``r n_cur_dam} \right)
= {1\over 2}\left( r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam] \right)
= r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] \notag
\end{equation}
n_cur_ani <- 8
n_cur_sire <- tbl_ped_exp8p1ext$Sire[n_cur_ani]
n_cur_dam <- tbl_ped_exp8p1ext$Dam[n_cur_ani]
matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] <- 0.5 * (matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam])
\begin{equation}
(A){r n_ex8p1_first_animal``r n_cur_ani} = {1\over 2}\left( (A){r n_ex8p1_first_animal``r n_cur_sire} + (A)_{r n_ex8p1_first_animal``r n_cur_dam} \right)
= {1\over 2}\left( r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam] \right)
= r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] \notag
\end{equation}
n_cur_ani <- 9
n_cur_sire <- tbl_ped_exp8p1ext$Sire[n_cur_ani]
n_cur_dam <- tbl_ped_exp8p1ext$Dam[n_cur_ani]
matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] <- 0.5 * (matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam])
\begin{equation}
(A){r n_ex8p1_first_animal``r n_cur_ani} = {1\over 2}\left( (A){r n_ex8p1_first_animal``r n_cur_sire} + (A)_{r n_ex8p1_first_animal``r n_cur_dam} \right)
= {1\over 2}\left( r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam] \right)
= r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] \notag
\end{equation}
n_cur_ani <- 10
n_cur_sire <- tbl_ped_exp8p1ext$Sire[n_cur_ani]
n_cur_dam <- tbl_ped_exp8p1ext$Dam[n_cur_ani]
matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] <- 0.5 * (matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam])
\begin{equation}
(A){r n_ex8p1_first_animal``r n_cur_ani} = {1\over 2}\left( (A){r n_ex8p1_first_animal``r n_cur_sire} + (A)_{r n_ex8p1_first_animal``r n_cur_dam} \right)
= {1\over 2}\left( r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam] \right)
= r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] \notag
\end{equation}
As a result, we have the first row of $A$
matA_empty_string_ex8p1[n_ex8p1_first_animal,] <- matA_empty_ex8p1[n_ex8p1_first_animal,]
### # display
cat(paste0(rmdhelp::bmatrix(pmat = matA_empty_string_ex8p1, ps_name = 'A', ps_env = '$$')))
Step 5
Copy the first row to the first column
matA_empty_ex8p1[2:n_nr_ani_ex8p1pedext,1] <- matA_empty_ex8p1[1, 2:n_nr_ani_ex8p1pedext]
matA_empty_string_ex8p1[2:n_nr_ani_ex8p1pedext,1] <- matA_empty_ex8p1[2:n_nr_ani_ex8p1pedext,1]
### # display
cat(paste0(rmdhelp::bmatrix(pmat = matA_empty_string_ex8p1, ps_name = 'A', ps_env = '$$')))
Step 6
Continue the same way with rows $2$ to $10$
Verification
We first have to specify the pedigree, before being able to get the numerator relationship matrix
n_nr_animals <- 10
suppressPackageStartupMessages( library(pedigreemm) )
ped <- pedigree(sire = c(NA,NA,NA,NA,1,1,4,4,4,4),
dam = c(NA,NA,NA,NA,2,3,5,5,6,6),
label = as.character(1:n_nr_animals))
mata_ex8p1_verify <- getA(ped = ped)
mata_ex8p1_verify
In the above result all elements which are $0$ are represented by a dot.
Problem 2: BLUP Animal Model
Use the following dataset to predict breeding values for all animals.
n_sigmae2 <- 24
n_sigmaa2 <- 8
s_ex8p2_data <- "https://charlotte-ngs.github.io/lbgfs2020/misc/data_ex06_p02.csv"
tbl_ex8p2_data <- readr::read_csv(file = s_ex8p2_data)
knitr::kable(tbl_ex8p2_data,
booktabs = TRUE,
longtable = TRUE,
caption = "Data for Animal Model")
Assumptions
- Random residuals are un-correlated and they all have equal variance $\sigma_e^2$ which is assumed to be
r n_sigmae2.
- The additive genetic variance $\sigma_u^2$ is assumed to be
r n_sigmaa2.
- The pedigree is the same as in Problem 1. You can use
solve() in R or pedigreemm::getAInv() to invert $A$.
Your Tasks
- Specify all components including expected values and variances of the animal model using the information from the dataset.
- Set up mixed model equations
- Solve mixed model equations for estimates of fixed effects and for predicted breeding values
Solution
The animal model in general has the following form
$$y = X \beta + Zu + e$$
\begin{tabular}{lll}
where & & \
& $y$ & vector of length $n$ of observations \
& $\beta$ & vector of length $p$ of unknown fixed effects \
& $X$ & $n\times p$ design matrix linking fixed effects to observations \
& $u$ & vector of length $q$ of unkown random breeding values \
& $Z$ & $n\times q$ design matrix linking breeding values to observations \
& $e$ & vector of length $n$ of unknown random residuals \
\end{tabular}
The expected values of the fixed effects $\beta$ are the fixed effects themselves, hence $E(\beta) = \beta$. The expected values of the random components are defined as
\begin{align}
E(u) &= 0 \notag \
E(e) &= 0 \notag \
E(y) &= X\beta \notag
\end{align}
The variances of fixed effects are always $0$. Based on the assumption of uncorrelated residuals, we know that $var(e) = I\sigma_e^2$. Because, we have an animal model, we also know that $var(u) = A \sigma_u^2$ where $A$ corresponds to the numerator relationship matrix. In summary, the variances of the random effects are
\begin{align}
var(u) &= G = A \sigma_u^2 \notag \
var(e) &= R = I \sigma_e^2 \notag \
var(y) &= ZGZ^T + R \notag
\end{align}
Inserting the information from the dataset into the model gives the following results.
- Vector of observations
vec_y <- as.vector(tbl_ex8p2_data$Observation)
cat(paste(rmdhelp::bcolumn_vector(pvec = vec_y, ps_name = 'y', ps_env = '$$')))
- Herds as fixed effects. We have two herds, hence vector $\beta$ has length $p = 2$. Component $\beta_1$ will denote the effect of the first herd and component $\beta_2$ the effect of the second herd.
n_nr_obs <- length(vec_y)
mat_x <- matrix(c(1, 0,
1, 0,
1, 0,
0, 1,
0, 1,
0, 1), nrow = n_nr_obs, byrow = TRUE)
vec_beta <- rmdhelp::vecGetVecElem(psBaseElement = "\\beta", pnVecLen = ncol(mat_x), psResult = "latex")
cat("$$\n")
cat(paste(rmdhelp::bmatrix(pmat = mat_x, ps_name = 'X')))
cat(" \\text{, } \n")
cat(paste(rmdhelp::bcolumn_vector(pvec = vec_beta, ps_name = '\\beta')))
cat("$$\n")
n_nr_animals <- 10
- Breeding values as random effects. In total, there are
r n_nr_animals animals in the pedigree and hence the length of the vector $a$ is $q = r n_nr_animals$.
mat_z <- cbind(matrix(0, nrow = n_nr_obs, ncol = n_nr_founder_ex8p1pedext), diag(1, nrow = n_nr_obs, ncol = n_nr_obs))
vec_a <- rmdhelp::vecGetVecElem(psBaseElement = "a", pnVecLen = n_nr_animals, psResult = "latex")
cat("$$\n")
cat(paste(rmdhelp::bmatrix(pmat = mat_z, ps_name = 'Z')))
cat(" \\text{, } \n")
cat(paste(rmdhelp::bcolumn_vector(pvec = vec_a, ps_name = 'u')))
cat("$$\n")
- The vector of random residuals is just
vec_e <- rmdhelp::vecGetVecElem(psBaseElement = "e", pnVecLen = n_nr_obs, psResult = "latex")
cat(paste(rmdhelp::bcolumn_vector(pvec = vec_e, ps_name = 'e', ps_env = '$$')))
The solutions for $\hat{\beta}$ and $\hat{a}$ are obtained by solving the mixed model equations. The mixed model equations for the animal model and under the assumptions specified above are defined as
$$
\left[
\begin{array}{lr}
X^TX & X^TZ \
Z^TX & Z^TZ + \lambda * A^{-1}
\end{array}
\right]
\left[
\begin{array}{c}
\hat{\beta} \
\hat{a}
\end{array}
\right]
=
\left[
\begin{array}{c}
X^Ty \
Z^Ty
\end{array}
\right]
$$
The single components are computed as
### # design matrices
mat_xtx <- crossprod(mat_x)
mat_xtz <- crossprod(mat_x, mat_z)
mat_ztz <- crossprod(mat_z)
### # variance ratio
lambda <- n_sigmae2 / n_sigmaa2
### # pedigree
suppressPackageStartupMessages( library(pedigreemm) )
ped <- pedigree(sire = c(rep(NA, (n_nr_animals-n_nr_obs)), tbl_ex8p2_data$Sire),
dam = c(rep(NA, (n_nr_animals-n_nr_obs)), tbl_ex8p2_data$Dam),
label = as.character(1:n_nr_animals))
mat_ainv = as.matrix(getAInv(ped = ped))
mat_ztzlainv = mat_ztz + lambda * mat_ainv
### # coefficient matrix
mat_coef <- rbind(cbind(mat_xtx, mat_xtz),cbind(t(mat_xtz), mat_ztzlainv))
### # unknows
vec_unknowns <- c("\\hat{\\beta}", "\\hat{a}")
### # right hand side
mat_xty <- crossprod(mat_x, vec_y)
mat_zty <- crossprod(mat_z, vec_y)
mat_rhs <- rbind(mat_xty, mat_zty)
### # solution
mat_sol <- solve(mat_coef, mat_rhs)
### # show the components
### # X^TX
cat("$$\n")
cat(paste(rmdhelp::bmatrix(pmat = mat_xtx, ps_name = 'X^TX')))
cat(" \\text{, } \n")
### # X^TZ
cat(paste(rmdhelp::bmatrix(pmat = mat_xtz, ps_name = 'X^TZ')))
cat("$$\n\n")
### # Z^TZ
cat(paste(rmdhelp::bmatrix(pmat = mat_ztz, ps_name = 'Z^TZ', ps_env = '$$')))
### # A^{-1}
cat(paste(rmdhelp::bmatrix(pmat = mat_ainv, ps_name = 'A^{-1}', ps_env = '$$')))
### # right hand side
cat("$$\n")
cat(paste(rmdhelp::bmatrix(pmat = mat_xty, ps_name = 'X^Ty')))
cat(" \\text{, } \n")
cat(paste(rmdhelp::bmatrix(pmat = mat_zty, ps_name = 'Z^Ty')))
cat("$$\n")
Putting everything together into the mixed model equations leads to the following results
cat("$$\n")
# cat(paste(rmdhelp::bmatrix(pmat = mat_coef)), "\n")
cat("\\cdot")
cat(paste(rmdhelp::bcolumn_vector(pvec = vec_unknowns)), '\n')
cat(" = \n")
cat(paste(rmdhelp::bmatrix(pmat = mat_rhs)))
cat("$$\n")
The solutions are computed as
cat("$$\n")
cat(paste(rmdhelp::bcolumn_vector(pvec = vec_unknowns)))
cat(" = \n")
cat(paste(rmdhelp::bmatrix(pmat = mat_coef)))
cat("^{-1}\n")
cat(paste(rmdhelp::bmatrix(pmat = mat_rhs)))
cat("$$\n")
The solutions are
vec_beta_hat <- rmdhelp::vecGetVecElem(psBaseElement = "\\hat{\\beta}", pnVecLen = ncol(mat_x), psResult = "latex")
vec_a_hat <- rmdhelp::vecGetVecElem(psBaseElement = "\\hat{a}", pnVecLen = n_nr_animals, psResult = "latex")
cat("$$\n")
cat(paste(rmdhelp::bcolumn_vector(pvec = c(vec_beta_hat, vec_a_hat))), '\n')
cat(" = \n")
cat(paste(rmdhelp::bmatrix(pmat = mat_sol)))
cat("$$\n")
charlotte-ngs/lbgfs2020 documentation built on Dec. 20, 2020, 5:39 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = FALSE, results = 'asis', fig.pos = 'H') knitr::knit_hooks$set(hook_convert_odg = rmdhelp::hook_convert_odg)
Problem 1: Numerator Relationship Matrix
Construct the numerator relationship matrix $A$ for the following pedigree and verify the result using the function getA() from package pedigreemm.
tbl_ped_ex8p1 <- tibble::tibble(Animal = c(5, 6, 7, 8, 9, 10), Sire = c(1, 1, 4, 4, 4, 4), Dam = c(2, 3, 5, 5, 6, 6)) knitr::kable(tbl_ped_ex8p1, booktabs = TRUE, longtable = TRUE, caption = "Pedigree For Constructing Numerator Relationship Matrix")
Solution
The numerator relationship is constructed using the follownig step-wise procedure. The following rules are used to compute the single elements.
-
Case 1: If both parents $s$ and $d$ of animal $i$ are known then
- the diagonal element $(A){ii}$ corresponds to: $(A){ii} = 1 + F_i = 1 + {1\over 2} (A)_{sd}$ and
- the offdiagonal element $(A){ji}$ is computed as: $(A){ji} = {1\over 2} ((A){js} + (A){jd})$
- because $A$ is symmetric $(A){ji} = (A){ij}$
-
Case 2: If only one parent $s$ is known and assumed unrelated to the mate
- $(A)_{ii} = 1$
- $(A){ij} = (A){ji} = {1\over 2} ((A)_{js}$
-
Case 3: If both parents are unknown
- $(A)_{ii} = 1$
- $(A){ij} = (A){ji} = 0$
Step 1
First, we extend the pedigree given in Table \@ref(tab:ex8p1ped). All animals without parents are added at the top of the pedigree. This results in the matrix shown in Table \@ref(tab:ex8p1pedext).
n_nr_ani_ex8p1pedext <- max(tbl_ped_ex8p1$Animal) n_nr_founder_ex8p1pedext <- min(tbl_ped_ex8p1$Animal) - 1 tbl_ped_exp8p1ext <- tibble::tibble(Animal = c(1:n_nr_ani_ex8p1pedext), Sire = c(rep(NA, n_nr_founder_ex8p1pedext), tbl_ped_ex8p1$Sire), Dam = c(rep(NA, n_nr_founder_ex8p1pedext), tbl_ped_ex8p1$Dam)) knitr::kable(tbl_ped_exp8p1ext, booktabs = TRUE, longtable = TRUE, caption = "Extended Pedigree")
Because the pedigree in Table \@ref(tab:ex8p1pedext) is already ordered such that parents are before offspring, we can directly go to the next step.
Step 2
We start with an empty numerator relationship matrix $A$. The matrix $A$ has dimensions $r n_nr_ani_ex8p1pedext \times r n_nr_ani_ex8p1pedext$
matA_empty_ex8p1 <- matrix(NA, nrow = n_nr_ani_ex8p1pedext, ncol = n_nr_ani_ex8p1pedext) matA_empty_string_ex8p1 <- matrix('', nrow = n_nr_ani_ex8p1pedext, ncol = n_nr_ani_ex8p1pedext) ### # display cat(paste(rmdhelp::bmatrix(pmat = matA_empty_string_ex8p1, ps_name = 'A', ps_env = '$$')))
Step 3
The single elements of $A$ are computed according to the rules listed above.
n_ex8p1_first_animal <- tbl_ped_exp8p1ext$Animal[1]
The computation is started with animal $r n_ex8p1_first_animal$. The first element is always the diagonal-element that corresponds to animal that we are currently looking at. For animal $r n_ex8p1_first_animal$ the diagonal element is $(A)_{r n_ex8p1_first_animal``r n_ex8p1_first_animal}$. Because animal r n_ex8p1_first_animal has not parents, we are in case 3 for the diagonal element. If an animal has unknown parents, it also means that the animals's inbreeding coefficient $F_i$ is $0$. Hence
$$(A)_{r n_ex8p1_first_animal``r n_ex8p1_first_animal} = 1$$
Now we have the first element of our numerator relationship matrix.
matA_empty_ex8p1[n_ex8p1_first_animal,n_ex8p1_first_animal] <- 1 matA_empty_string_ex8p1[n_ex8p1_first_animal,n_ex8p1_first_animal] <- 1 ### # display cat(paste0(rmdhelp::bmatrix(pmat = matA_empty_string_ex8p1, ps_name = 'A', ps_env = '$$')))
Step 4
The next elements that need to be computed are the off-diagnoal element on row $1$. Elements $A_{12}$, $A_{13}$ and $A_{14}$ correspond to the additive genetic relationship between animal $1$ and animals $2$, $3$ and $4$. Because animals $2$, $3$ and $4$ all have unknown parents, we are for all three elements in case 3, hence we can state
$$(A){12} = (A){13} = (A)_{14} = 0$$
matA_empty_ex8p1[n_ex8p1_first_animal,2:4] <- 0
For the remaining elements of the first row of $A$, the elements correspond to the additive genetic relationship between animal $1$ and animals $5$ to $10$. Because animals $5$ to $10$ all have known parents, we have to use case 1 in the above formulated rules.
n_cur_ani <- 5 n_cur_sire <- tbl_ped_exp8p1ext$Sire[n_cur_ani] n_cur_dam <- tbl_ped_exp8p1ext$Dam[n_cur_ani] matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] <- 0.5 * (matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam])
\begin{equation}
(A){r n_ex8p1_first_animal``r n_cur_ani} = {1\over 2}\left( (A){r n_ex8p1_first_animal``r n_cur_sire} + (A)_{r n_ex8p1_first_animal``r n_cur_dam} \right)
= {1\over 2}\left( r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam] \right)
= r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] \notag
\end{equation}
n_cur_ani <- 6 n_cur_sire <- tbl_ped_exp8p1ext$Sire[n_cur_ani] n_cur_dam <- tbl_ped_exp8p1ext$Dam[n_cur_ani] matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] <- 0.5 * (matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam])
\begin{equation}
(A){r n_ex8p1_first_animal``r n_cur_ani} = {1\over 2}\left( (A){r n_ex8p1_first_animal``r n_cur_sire} + (A)_{r n_ex8p1_first_animal``r n_cur_dam} \right)
= {1\over 2}\left( r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam] \right)
= r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] \notag
\end{equation}
n_cur_ani <- 7 n_cur_sire <- tbl_ped_exp8p1ext$Sire[n_cur_ani] n_cur_dam <- tbl_ped_exp8p1ext$Dam[n_cur_ani] matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] <- 0.5 * (matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam])
\begin{equation}
(A){r n_ex8p1_first_animal``r n_cur_ani} = {1\over 2}\left( (A){r n_ex8p1_first_animal``r n_cur_sire} + (A)_{r n_ex8p1_first_animal``r n_cur_dam} \right)
= {1\over 2}\left( r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam] \right)
= r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] \notag
\end{equation}
n_cur_ani <- 8 n_cur_sire <- tbl_ped_exp8p1ext$Sire[n_cur_ani] n_cur_dam <- tbl_ped_exp8p1ext$Dam[n_cur_ani] matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] <- 0.5 * (matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam])
\begin{equation}
(A){r n_ex8p1_first_animal``r n_cur_ani} = {1\over 2}\left( (A){r n_ex8p1_first_animal``r n_cur_sire} + (A)_{r n_ex8p1_first_animal``r n_cur_dam} \right)
= {1\over 2}\left( r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam] \right)
= r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] \notag
\end{equation}
n_cur_ani <- 9 n_cur_sire <- tbl_ped_exp8p1ext$Sire[n_cur_ani] n_cur_dam <- tbl_ped_exp8p1ext$Dam[n_cur_ani] matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] <- 0.5 * (matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam])
\begin{equation}
(A){r n_ex8p1_first_animal``r n_cur_ani} = {1\over 2}\left( (A){r n_ex8p1_first_animal``r n_cur_sire} + (A)_{r n_ex8p1_first_animal``r n_cur_dam} \right)
= {1\over 2}\left( r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam] \right)
= r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] \notag
\end{equation}
n_cur_ani <- 10 n_cur_sire <- tbl_ped_exp8p1ext$Sire[n_cur_ani] n_cur_dam <- tbl_ped_exp8p1ext$Dam[n_cur_ani] matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] <- 0.5 * (matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam])
\begin{equation}
(A){r n_ex8p1_first_animal``r n_cur_ani} = {1\over 2}\left( (A){r n_ex8p1_first_animal``r n_cur_sire} + (A)_{r n_ex8p1_first_animal``r n_cur_dam} \right)
= {1\over 2}\left( r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_sire] + r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_dam] \right)
= r matA_empty_ex8p1[n_ex8p1_first_animal,n_cur_ani] \notag
\end{equation}
As a result, we have the first row of $A$
matA_empty_string_ex8p1[n_ex8p1_first_animal,] <- matA_empty_ex8p1[n_ex8p1_first_animal,] ### # display cat(paste0(rmdhelp::bmatrix(pmat = matA_empty_string_ex8p1, ps_name = 'A', ps_env = '$$')))
Step 5
Copy the first row to the first column
matA_empty_ex8p1[2:n_nr_ani_ex8p1pedext,1] <- matA_empty_ex8p1[1, 2:n_nr_ani_ex8p1pedext] matA_empty_string_ex8p1[2:n_nr_ani_ex8p1pedext,1] <- matA_empty_ex8p1[2:n_nr_ani_ex8p1pedext,1] ### # display cat(paste0(rmdhelp::bmatrix(pmat = matA_empty_string_ex8p1, ps_name = 'A', ps_env = '$$')))
Step 6
Continue the same way with rows $2$ to $10$
Verification
We first have to specify the pedigree, before being able to get the numerator relationship matrix
n_nr_animals <- 10 suppressPackageStartupMessages( library(pedigreemm) ) ped <- pedigree(sire = c(NA,NA,NA,NA,1,1,4,4,4,4), dam = c(NA,NA,NA,NA,2,3,5,5,6,6), label = as.character(1:n_nr_animals)) mata_ex8p1_verify <- getA(ped = ped) mata_ex8p1_verify
In the above result all elements which are $0$ are represented by a dot.
Problem 2: BLUP Animal Model
Use the following dataset to predict breeding values for all animals.
n_sigmae2 <- 24 n_sigmaa2 <- 8 s_ex8p2_data <- "https://charlotte-ngs.github.io/lbgfs2020/misc/data_ex06_p02.csv" tbl_ex8p2_data <- readr::read_csv(file = s_ex8p2_data) knitr::kable(tbl_ex8p2_data, booktabs = TRUE, longtable = TRUE, caption = "Data for Animal Model")
Assumptions
- Random residuals are un-correlated and they all have equal variance $\sigma_e^2$ which is assumed to be
r n_sigmae2. - The additive genetic variance $\sigma_u^2$ is assumed to be
r n_sigmaa2. - The pedigree is the same as in Problem 1. You can use
solve()in R orpedigreemm::getAInv()to invert $A$.
Your Tasks
- Specify all components including expected values and variances of the animal model using the information from the dataset.
- Set up mixed model equations
- Solve mixed model equations for estimates of fixed effects and for predicted breeding values
Solution
The animal model in general has the following form
$$y = X \beta + Zu + e$$
\begin{tabular}{lll} where & & \ & $y$ & vector of length $n$ of observations \ & $\beta$ & vector of length $p$ of unknown fixed effects \ & $X$ & $n\times p$ design matrix linking fixed effects to observations \ & $u$ & vector of length $q$ of unkown random breeding values \ & $Z$ & $n\times q$ design matrix linking breeding values to observations \ & $e$ & vector of length $n$ of unknown random residuals \ \end{tabular}
The expected values of the fixed effects $\beta$ are the fixed effects themselves, hence $E(\beta) = \beta$. The expected values of the random components are defined as
\begin{align} E(u) &= 0 \notag \ E(e) &= 0 \notag \ E(y) &= X\beta \notag \end{align}
The variances of fixed effects are always $0$. Based on the assumption of uncorrelated residuals, we know that $var(e) = I\sigma_e^2$. Because, we have an animal model, we also know that $var(u) = A \sigma_u^2$ where $A$ corresponds to the numerator relationship matrix. In summary, the variances of the random effects are
\begin{align} var(u) &= G = A \sigma_u^2 \notag \ var(e) &= R = I \sigma_e^2 \notag \ var(y) &= ZGZ^T + R \notag \end{align}
Inserting the information from the dataset into the model gives the following results.
- Vector of observations
vec_y <- as.vector(tbl_ex8p2_data$Observation) cat(paste(rmdhelp::bcolumn_vector(pvec = vec_y, ps_name = 'y', ps_env = '$$')))
- Herds as fixed effects. We have two herds, hence vector $\beta$ has length $p = 2$. Component $\beta_1$ will denote the effect of the first herd and component $\beta_2$ the effect of the second herd.
n_nr_obs <- length(vec_y) mat_x <- matrix(c(1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1), nrow = n_nr_obs, byrow = TRUE) vec_beta <- rmdhelp::vecGetVecElem(psBaseElement = "\\beta", pnVecLen = ncol(mat_x), psResult = "latex") cat("$$\n") cat(paste(rmdhelp::bmatrix(pmat = mat_x, ps_name = 'X'))) cat(" \\text{, } \n") cat(paste(rmdhelp::bcolumn_vector(pvec = vec_beta, ps_name = '\\beta'))) cat("$$\n")
n_nr_animals <- 10
- Breeding values as random effects. In total, there are
r n_nr_animalsanimals in the pedigree and hence the length of the vector $a$ is $q =r n_nr_animals$.
mat_z <- cbind(matrix(0, nrow = n_nr_obs, ncol = n_nr_founder_ex8p1pedext), diag(1, nrow = n_nr_obs, ncol = n_nr_obs)) vec_a <- rmdhelp::vecGetVecElem(psBaseElement = "a", pnVecLen = n_nr_animals, psResult = "latex") cat("$$\n") cat(paste(rmdhelp::bmatrix(pmat = mat_z, ps_name = 'Z'))) cat(" \\text{, } \n") cat(paste(rmdhelp::bcolumn_vector(pvec = vec_a, ps_name = 'u'))) cat("$$\n")
- The vector of random residuals is just
vec_e <- rmdhelp::vecGetVecElem(psBaseElement = "e", pnVecLen = n_nr_obs, psResult = "latex") cat(paste(rmdhelp::bcolumn_vector(pvec = vec_e, ps_name = 'e', ps_env = '$$')))
The solutions for $\hat{\beta}$ and $\hat{a}$ are obtained by solving the mixed model equations. The mixed model equations for the animal model and under the assumptions specified above are defined as
$$ \left[ \begin{array}{lr} X^TX & X^TZ \ Z^TX & Z^TZ + \lambda * A^{-1} \end{array} \right] \left[ \begin{array}{c} \hat{\beta} \ \hat{a} \end{array} \right] = \left[ \begin{array}{c} X^Ty \ Z^Ty \end{array} \right] $$
The single components are computed as
### # design matrices mat_xtx <- crossprod(mat_x) mat_xtz <- crossprod(mat_x, mat_z) mat_ztz <- crossprod(mat_z) ### # variance ratio lambda <- n_sigmae2 / n_sigmaa2 ### # pedigree suppressPackageStartupMessages( library(pedigreemm) ) ped <- pedigree(sire = c(rep(NA, (n_nr_animals-n_nr_obs)), tbl_ex8p2_data$Sire), dam = c(rep(NA, (n_nr_animals-n_nr_obs)), tbl_ex8p2_data$Dam), label = as.character(1:n_nr_animals)) mat_ainv = as.matrix(getAInv(ped = ped)) mat_ztzlainv = mat_ztz + lambda * mat_ainv ### # coefficient matrix mat_coef <- rbind(cbind(mat_xtx, mat_xtz),cbind(t(mat_xtz), mat_ztzlainv)) ### # unknows vec_unknowns <- c("\\hat{\\beta}", "\\hat{a}") ### # right hand side mat_xty <- crossprod(mat_x, vec_y) mat_zty <- crossprod(mat_z, vec_y) mat_rhs <- rbind(mat_xty, mat_zty) ### # solution mat_sol <- solve(mat_coef, mat_rhs) ### # show the components ### # X^TX cat("$$\n") cat(paste(rmdhelp::bmatrix(pmat = mat_xtx, ps_name = 'X^TX'))) cat(" \\text{, } \n") ### # X^TZ cat(paste(rmdhelp::bmatrix(pmat = mat_xtz, ps_name = 'X^TZ'))) cat("$$\n\n") ### # Z^TZ cat(paste(rmdhelp::bmatrix(pmat = mat_ztz, ps_name = 'Z^TZ', ps_env = '$$'))) ### # A^{-1} cat(paste(rmdhelp::bmatrix(pmat = mat_ainv, ps_name = 'A^{-1}', ps_env = '$$'))) ### # right hand side cat("$$\n") cat(paste(rmdhelp::bmatrix(pmat = mat_xty, ps_name = 'X^Ty'))) cat(" \\text{, } \n") cat(paste(rmdhelp::bmatrix(pmat = mat_zty, ps_name = 'Z^Ty'))) cat("$$\n")
Putting everything together into the mixed model equations leads to the following results
cat("$$\n") # cat(paste(rmdhelp::bmatrix(pmat = mat_coef)), "\n") cat("\\cdot") cat(paste(rmdhelp::bcolumn_vector(pvec = vec_unknowns)), '\n') cat(" = \n") cat(paste(rmdhelp::bmatrix(pmat = mat_rhs))) cat("$$\n")
The solutions are computed as
cat("$$\n") cat(paste(rmdhelp::bcolumn_vector(pvec = vec_unknowns))) cat(" = \n") cat(paste(rmdhelp::bmatrix(pmat = mat_coef))) cat("^{-1}\n") cat(paste(rmdhelp::bmatrix(pmat = mat_rhs))) cat("$$\n")
The solutions are
vec_beta_hat <- rmdhelp::vecGetVecElem(psBaseElement = "\\hat{\\beta}", pnVecLen = ncol(mat_x), psResult = "latex") vec_a_hat <- rmdhelp::vecGetVecElem(psBaseElement = "\\hat{a}", pnVecLen = n_nr_animals, psResult = "latex") cat("$$\n") cat(paste(rmdhelp::bcolumn_vector(pvec = c(vec_beta_hat, vec_a_hat))), '\n') cat(" = \n") cat(paste(rmdhelp::bmatrix(pmat = mat_sol))) cat("$$\n")
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.