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)
cnt <- rmdhelp::R6ClassCount$new()
cnt$set_prefix(ps_prefix = "## Problem")
cat(cnt$out(ps_suffix = "Multivariate BLUP Animal Model"), "\n")
n_nr_trait <- 2
n_nr_founder <- 3
n_nr_animal <- 8
n_nr_observation <- n_nr_animal - n_nr_founder
tbl_data_sol12p01 <- tibble::tibble(Animal = c((n_nr_founder+1):n_nr_animal),
Sex = c("Male", "Female","Female","Male","Male"),
Sire = c(1,3,1,4,3),
Dam = c(NA,2,2,5,6),
WWG = c(4.5,2.9,3.9,3.5,5.0),
PWG = c(6.8,5.0,6.8,6.0,7.5))
The table below contains data for pre-weaning gain (WWG) and post-weaning gain (PWG) for r n_nr_observation beef calves.
knitr::kable(tbl_data_sol12p01,
booktabs = TRUE,
longtable = TRUE)
The genetic variance-covariance matrix $G_0$ between the traits is
mat_g0 <- matrix(data = c(20,18,18,40), nrow = n_nr_trait, byrow = TRUE)
cat(paste(rmdhelp::bmatrix(pmat = mat_g0, ps_name = 'G_0', ps_env = '$$'), collapse = '\n'), '\n')
The residual variance-covariance matrix $R_0$ between the traits is
mat_r0 <- matrix(data = c(40,11,11,30), nrow = n_nr_trait, byrow = TRUE)
cat(paste(rmdhelp::bmatrix(pmat = mat_r0, ps_name = 'R_0', ps_env = '$$'), collapse = '\n'), '\n')
Your Task
Set up the mixed model equations for a multivariate BLUP analysis and compute the estimates for the fixed effects and the predictions for the breeding values.
Solution
The matrices $X_1$ and $X_2$ relate records of PWG and WWG to sex effects. For both traits, we have an effect for the male and female sex. Hence the vector $\beta$ of fixed effects corresponds to
vec_beta <- c("\\beta_{M,WWG}", "\\beta_{F,WWG}", "\\beta_{M,PWG}", "\\beta_{F,PWG}")
cat(paste(rmdhelp::bcolumn_vector(pvec = vec_beta, ps_name = '\\beta', ps_env = '$$'), collapse = '\n'), '\n')
The matrices $X_1$ and $X_2$ are the same and correspond to
mat_x1 <- mat_x2 <- matrix(data = c(1, 0,
0, 1,
0, 1,
1, 0,
1, 0), nrow = n_nr_observation, byrow = TRUE)
cat(paste(rmdhelp::bmatrix(pmat = mat_x1, ps_name = 'X_1 = X_2', ps_env = '$$'), collapse = '\n'), '\n')
Combining them to the multivariate version leads to
$$X = \left[
\begin{array}{lr}
X_1 & 0 \
0 & X_2
\end{array}
\right]$$
mat_zero <- matrix(0, nrow = nrow(mat_x1), ncol = ncol(mat_x1))
mat_x <- rbind(cbind(mat_x1, mat_zero), cbind(mat_zero, mat_x2))
cat(paste(rmdhelp::bmatrix(pmat = mat_x, ps_name = 'X', ps_env = '$$'), collapse = '\n'), '\n')
Using the matrix $X$ together with matrix $R = I_n \otimes R_0$ to get
mat_r <- mat_r0 %x% diag(1, n_nr_observation)
mat_rinv <- solve(mat_r)
mat_xtrinvx <- t(mat_x) %*% mat_rinv %*% mat_x
cat(paste(rmdhelp::bmatrix(pmat = mat_xtrinvx, ps_name = 'X^TR^{-1}X', ps_env = '$$'), collapse = '\n'), '\n')
Similarly to the fixed effects, we can put together the vector of breeding values $a$.
vec_a1 <- c("u_{1,WWG}",
"u_{2,WWG}",
"u_{3,WWG}",
"u_{4,WWG}",
"u_{5,WWG}",
"u_{6,WWG}",
"u_{7,WWG}",
"u_{8,WWG}")
vec_a2 <- c("u_{1,PWG}",
"u_{2,PWG}",
"u_{3,PWG}",
"u_{4,PWG}",
"u_{5,PWG}",
"u_{6,PWG}",
"u_{7,PWG}",
"u_{8,PWG}")
vec_a <- c(vec_a1, vec_a2)
cat(paste(rmdhelp::bcolumn_vector(pvec = vec_a, ps_name = 'u', ps_env = '$$'), collapse = '\n'), '\n')
The design matrices $Z_1$ and $Z_2$ are equal and they link observations to breeding values.
mat_z1zero <- matrix(0, nrow = n_nr_observation, ncol = n_nr_founder)
mat_z1 <- mat_z2 <- cbind(mat_z1zero, diag(1, n_nr_observation))
cat(paste(rmdhelp::bmatrix(pmat = mat_z1, ps_name = 'Z_1 = Z_2', ps_env = '$$'), collapse = '\n'), '\n')
$$Z = \left[
\begin{array}{lr}
Z_1 & 0 \
0 & Z_2
\end{array}
\right]$$
mat_zzero <- matrix(0, nrow = nrow(mat_z1), ncol(mat_z2))
mat_z <- rbind(cbind(mat_z1, mat_zzero), cbind(mat_zzero, mat_z2))
### # TODO: The following must be integrated into rmdhelp
### # ==> see https://tex.stackexchange.com/questions/95162/how-to-create-a-matrix-with-20-columns-in-latex/95163
cat('\\setcounter{MaxMatrixCols}{30}\n')
cat(paste(rmdhelp::bmatrix(pmat = mat_z, ps_name = 'Z', ps_env = '$$'), collapse = '\n'), '\n')
Together with the numerator relationship matrix $A$ we can get $G = G_0 \otimes A$ and from this $G^{-1} = G_0^{-1} \otimes A^{-1}$
ped_sol12p01 <- pedigreemm::pedigree(sire = c(rep(NA, n_nr_founder), tbl_data_sol12p01$Sire),
dam = c(rep(NA, n_nr_founder), tbl_data_sol12p01$Dam),
label = as.character(1:n_nr_animal))
mat_ainv <- as.matrix(pedigreemm::getAInv(ped = ped_sol12p01))
cat(paste(rmdhelp::bmatrix(pmat = round(mat_ainv, digits = 4), ps_name = 'A^{-1}', ps_env = '$$'), collapse = '\n'), '\n')
mat_ginv <- solve(mat_g0) %x% mat_ainv
cat("\\tiny \n")
cat(paste(rmdhelp::bmatrix(pmat = round(mat_ginv, digits = 3), ps_name = 'G^{-1}', ps_env = '$$'), collapse = '\n'), '\n')
cat("\\normalsize \n")
Using the matrics $X$, $Z$, $R^{-1}$ and $G^{-1}$, we can compute $Z^TR^{-1}X$ and $Z^TR^{-1}Z + G^{-1}$. These matrices define the coefficient matrix of the mixed model equations. But they are too be to be shown here.
The vector $y$ of observations contains all observations of both traits
vec_y1 <- tbl_data_sol12p01$WWG
vec_y2 <- tbl_data_sol12p01$PWG
vec_y <- c(vec_y1, vec_y2)
cat(paste(rmdhelp::bcolumn_vector(pvec = vec_y, ps_name = 'y', ps_env = '$$'), collapse = '\n'), '\n')
The right-hand side is computed as
$$
\left[
\begin{array}{c}
X^TR^{-1}y \
Z^TR^{-1}y
\end{array}
\right]
$$
The solutions are
### # coefficient matrix
mat_xtrinvz <- t(mat_x) %*% mat_rinv %*% mat_z
mat_ztrinvzginv <- t(mat_z) %*% mat_rinv %*% mat_z + mat_ginv
mat_coeff <- rbind(cbind(mat_xtrinvx, mat_xtrinvz), cbind(t(mat_xtrinvz), mat_ztrinvzginv))
### # right-hand side
mat_rhs <- rbind(t(mat_x) %*% mat_rinv %*% vec_y,
t(mat_z) %*% mat_rinv %*% vec_y)
vec_sol <- solve(mat_coeff, mat_rhs)
vec_beta_hat <- c("\\widehat{\\beta_{M,WWG}}",
"\\widehat{\\beta_{F,WWG}}",
"\\widehat{\\beta_{M,PWG}}",
"\\widehat{\\beta_{F,PWG}}")
vec_a1_hat <- c("\\widehat{u_{1,WWG}}",
"\\widehat{u_{2,WWG}}",
"\\widehat{u_{3,WWG}}",
"\\widehat{u_{4,WWG}}",
"\\widehat{u_{5,WWG}}",
"\\widehat{u_{6,WWG}}",
"\\widehat{u_{7,WWG}}",
"\\widehat{u_{8,WWG}}")
vec_a2_hat <- c("\\widehat{u_{1,PWG}}",
"\\widehat{u_{2,PWG}}",
"\\widehat{u_{3,PWG}}",
"\\widehat{u_{4,PWG}}",
"\\widehat{u_{5,PWG}}",
"\\widehat{u_{6,PWG}}",
"\\widehat{u_{7,PWG}}",
"\\widehat{u_{8,PWG}}")
vec_a_hat <- c(vec_a1_hat, vec_a2_hat)
vec_hat_unknown <- c(vec_beta_hat, vec_a_hat)
cat("$$\n")
cat(paste(rmdhelp::bcolumn_vector(pvec = vec_hat_unknown), collapse = '\n'), '\n')
cat(" = ")
cat(paste(rmdhelp::bmatrix(pmat = round(vec_sol, digits = 4)), collapse = '\n'), '\n')
cat("$$\n")
\vspace{3ex}
cat(cnt$out(ps_suffix = "Comparison of Reliabilites"), "\n")
Compare the predicted breeding values and the reliabilites obtained as results of Problem 1 with results from two univariate analyses for the same traits are used in Problem 1. All parameters can be taken from Problem 1.
Solution
For a predicted breeding value $\hat{u}_i$, the reliability $B_i$ is computed as
$$B_i = r_{u,\hat{u}}^2 = 1 - \frac{PEV(\hat{u}i)}{var(u_i)} = 1 - \frac{C{ii}^{22}}{var(u_i)}$$
where $C_{ii}^{22}$ are obtained from the inverse coefficient matrix of the mixed model equations. Just as a reminder, we can write the mixed model equations (MME) as
$$M \cdot s = r$$
with the vectors $r$ and $s$ corresponding to the right-hand side and to the unknowns of the MME. Hence
$$r = \left[\begin{array}{c}X^TR^{-1}y \ Z^TR^{-1}y \end{array}\right]$$
and
$$s = \left[\begin{array}{c}\hat{\beta} \ \hat{u} \end{array}\right]$$
The matrix $C^{22}$ is taken from the inverse coefficient matrix.
$$
M^{-1}
=
\left[
\begin{array}{lr}
X^TR^{-1}X & X^TR^{-1}Z \
Z^TR^{-1}X & Z^TR^{-1}Z + G^{-1}
\end{array}
\right]^{-1}
=
\left[
\begin{array}{lr}
C^{11} & C^{12} \
C^{21} & C^{22}
\end{array}
\right]
$$
For the two univariate analyses, we get the solutions for the fixed effects and the breeding values and their reliabilities as follows
- WWG: estimates $s_{WWG}$ and reliabilites $B_{WWG}$
lambda1 <- mat_r0[1,1] / mat_g0[1,1]
matxtx1 <- crossprod(mat_x1)
matxtz1 <- crossprod(mat_x1,mat_z1)
matztx1 <- t(matxtz1)
matztzainv1 <- crossprod(mat_z1) + mat_ainv * lambda1
mat_coeff1 <- rbind(cbind(matxtx1, matxtz1), cbind(matztx1, matztzainv1))
mat_rhs1 <- rbind(crossprod(mat_x1, vec_y1), crossprod(mat_z1, vec_y1))
mat_sol1 <- solve(mat_coeff1, mat_rhs1)
cat(paste(rmdhelp::bmatrix(pmat = round(mat_sol1, digits = 4), ps_name = 's_{WWG}', ps_env = '$$'), collapse = '\n'), '\n')
mat_coeff_inv1 <- solve(mat_coeff1 / mat_r0[1,1])
n_nr_sol <- dim(mat_coeff_inv1)[1]
vec_c22_diag1 <- diag(mat_coeff_inv1)[(n_nr_sol - n_nr_animal + 1):n_nr_sol]
vec_rel1 <- 1 - vec_c22_diag1 / mat_g0[1,1]
cat(paste(rmdhelp::bcolumn_vector(pvec = round(vec_rel1, digits = 4), ps_name = 'B_{WWG}', ps_env = '$$'), collapse = '\n'), '\n')
- PWG
lambda2 <- mat_r0[2,2] / mat_g0[2,2]
matxtx2 <- crossprod(mat_x2)
matxtz2 <- crossprod(mat_x2,mat_z2)
matztx2 <- t(matxtz2)
matztzainv2 <- crossprod(mat_z2) + mat_ainv * lambda2
mat_coeff2 <- rbind(cbind(matxtx2, matxtz2), cbind(matztx2, matztzainv2))
mat_rhs2 <- rbind(crossprod(mat_x2, vec_y2), crossprod(mat_z2, vec_y2))
mat_sol2 <- solve(mat_coeff2, mat_rhs2)
cat(paste(rmdhelp::bmatrix(pmat = round(mat_sol2, digits = 4), ps_name = 's_{PWG}', ps_env = '$$'), collapse = '\n'), '\n')
mat_coeff_inv2 <- solve(mat_coeff2 / mat_r0[2,2])
n_nr_sol <- dim(mat_coeff_inv2)[1]
vec_c22_diag2 <- diag(mat_coeff_inv2)[(n_nr_sol - n_nr_animal + 1):n_nr_sol]
vec_rel2 <- 1 - vec_c22_diag2 / mat_g0[2,2]
cat(paste(rmdhelp::bcolumn_vector(pvec = round(vec_rel2, digits = 4), ps_name = 'B_{PWG}', ps_env = '$$'), collapse = '\n'), '\n')
The reliabilities from the bivariate analysis are obtained as
mat_coeff_inv <- solve(mat_coeff)
n_nr_sol <- dim(mat_coeff_inv)[1]
vec_c22_diag <- diag(mat_coeff_inv)[(n_nr_sol - n_nr_trait * n_nr_animal + 1):n_nr_sol]
vec_rel <- 1 - vec_c22_diag / c(rep(mat_g0[1,1], n_nr_animal), rep(mat_g0[2,2], n_nr_animal))
cat(paste(rmdhelp::bcolumn_vector(pvec = round(vec_rel, digits = 4), ps_name = 'B', ps_env = '$$'), collapse = '\n'), '\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)
cnt <- rmdhelp::R6ClassCount$new() cnt$set_prefix(ps_prefix = "## Problem")
cat(cnt$out(ps_suffix = "Multivariate BLUP Animal Model"), "\n")
n_nr_trait <- 2 n_nr_founder <- 3 n_nr_animal <- 8 n_nr_observation <- n_nr_animal - n_nr_founder tbl_data_sol12p01 <- tibble::tibble(Animal = c((n_nr_founder+1):n_nr_animal), Sex = c("Male", "Female","Female","Male","Male"), Sire = c(1,3,1,4,3), Dam = c(NA,2,2,5,6), WWG = c(4.5,2.9,3.9,3.5,5.0), PWG = c(6.8,5.0,6.8,6.0,7.5))
The table below contains data for pre-weaning gain (WWG) and post-weaning gain (PWG) for r n_nr_observation beef calves.
knitr::kable(tbl_data_sol12p01, booktabs = TRUE, longtable = TRUE)
The genetic variance-covariance matrix $G_0$ between the traits is
mat_g0 <- matrix(data = c(20,18,18,40), nrow = n_nr_trait, byrow = TRUE) cat(paste(rmdhelp::bmatrix(pmat = mat_g0, ps_name = 'G_0', ps_env = '$$'), collapse = '\n'), '\n')
The residual variance-covariance matrix $R_0$ between the traits is
mat_r0 <- matrix(data = c(40,11,11,30), nrow = n_nr_trait, byrow = TRUE) cat(paste(rmdhelp::bmatrix(pmat = mat_r0, ps_name = 'R_0', ps_env = '$$'), collapse = '\n'), '\n')
Your Task
Set up the mixed model equations for a multivariate BLUP analysis and compute the estimates for the fixed effects and the predictions for the breeding values.
Solution
The matrices $X_1$ and $X_2$ relate records of PWG and WWG to sex effects. For both traits, we have an effect for the male and female sex. Hence the vector $\beta$ of fixed effects corresponds to
vec_beta <- c("\\beta_{M,WWG}", "\\beta_{F,WWG}", "\\beta_{M,PWG}", "\\beta_{F,PWG}") cat(paste(rmdhelp::bcolumn_vector(pvec = vec_beta, ps_name = '\\beta', ps_env = '$$'), collapse = '\n'), '\n')
The matrices $X_1$ and $X_2$ are the same and correspond to
mat_x1 <- mat_x2 <- matrix(data = c(1, 0, 0, 1, 0, 1, 1, 0, 1, 0), nrow = n_nr_observation, byrow = TRUE) cat(paste(rmdhelp::bmatrix(pmat = mat_x1, ps_name = 'X_1 = X_2', ps_env = '$$'), collapse = '\n'), '\n')
Combining them to the multivariate version leads to
$$X = \left[ \begin{array}{lr} X_1 & 0 \ 0 & X_2 \end{array} \right]$$
mat_zero <- matrix(0, nrow = nrow(mat_x1), ncol = ncol(mat_x1)) mat_x <- rbind(cbind(mat_x1, mat_zero), cbind(mat_zero, mat_x2)) cat(paste(rmdhelp::bmatrix(pmat = mat_x, ps_name = 'X', ps_env = '$$'), collapse = '\n'), '\n')
Using the matrix $X$ together with matrix $R = I_n \otimes R_0$ to get
mat_r <- mat_r0 %x% diag(1, n_nr_observation) mat_rinv <- solve(mat_r) mat_xtrinvx <- t(mat_x) %*% mat_rinv %*% mat_x cat(paste(rmdhelp::bmatrix(pmat = mat_xtrinvx, ps_name = 'X^TR^{-1}X', ps_env = '$$'), collapse = '\n'), '\n')
Similarly to the fixed effects, we can put together the vector of breeding values $a$.
vec_a1 <- c("u_{1,WWG}", "u_{2,WWG}", "u_{3,WWG}", "u_{4,WWG}", "u_{5,WWG}", "u_{6,WWG}", "u_{7,WWG}", "u_{8,WWG}") vec_a2 <- c("u_{1,PWG}", "u_{2,PWG}", "u_{3,PWG}", "u_{4,PWG}", "u_{5,PWG}", "u_{6,PWG}", "u_{7,PWG}", "u_{8,PWG}") vec_a <- c(vec_a1, vec_a2) cat(paste(rmdhelp::bcolumn_vector(pvec = vec_a, ps_name = 'u', ps_env = '$$'), collapse = '\n'), '\n')
The design matrices $Z_1$ and $Z_2$ are equal and they link observations to breeding values.
mat_z1zero <- matrix(0, nrow = n_nr_observation, ncol = n_nr_founder) mat_z1 <- mat_z2 <- cbind(mat_z1zero, diag(1, n_nr_observation)) cat(paste(rmdhelp::bmatrix(pmat = mat_z1, ps_name = 'Z_1 = Z_2', ps_env = '$$'), collapse = '\n'), '\n')
$$Z = \left[ \begin{array}{lr} Z_1 & 0 \ 0 & Z_2 \end{array} \right]$$
mat_zzero <- matrix(0, nrow = nrow(mat_z1), ncol(mat_z2)) mat_z <- rbind(cbind(mat_z1, mat_zzero), cbind(mat_zzero, mat_z2)) ### # TODO: The following must be integrated into rmdhelp ### # ==> see https://tex.stackexchange.com/questions/95162/how-to-create-a-matrix-with-20-columns-in-latex/95163 cat('\\setcounter{MaxMatrixCols}{30}\n') cat(paste(rmdhelp::bmatrix(pmat = mat_z, ps_name = 'Z', ps_env = '$$'), collapse = '\n'), '\n')
Together with the numerator relationship matrix $A$ we can get $G = G_0 \otimes A$ and from this $G^{-1} = G_0^{-1} \otimes A^{-1}$
ped_sol12p01 <- pedigreemm::pedigree(sire = c(rep(NA, n_nr_founder), tbl_data_sol12p01$Sire), dam = c(rep(NA, n_nr_founder), tbl_data_sol12p01$Dam), label = as.character(1:n_nr_animal)) mat_ainv <- as.matrix(pedigreemm::getAInv(ped = ped_sol12p01)) cat(paste(rmdhelp::bmatrix(pmat = round(mat_ainv, digits = 4), ps_name = 'A^{-1}', ps_env = '$$'), collapse = '\n'), '\n')
mat_ginv <- solve(mat_g0) %x% mat_ainv cat("\\tiny \n") cat(paste(rmdhelp::bmatrix(pmat = round(mat_ginv, digits = 3), ps_name = 'G^{-1}', ps_env = '$$'), collapse = '\n'), '\n') cat("\\normalsize \n")
Using the matrics $X$, $Z$, $R^{-1}$ and $G^{-1}$, we can compute $Z^TR^{-1}X$ and $Z^TR^{-1}Z + G^{-1}$. These matrices define the coefficient matrix of the mixed model equations. But they are too be to be shown here.
The vector $y$ of observations contains all observations of both traits
vec_y1 <- tbl_data_sol12p01$WWG vec_y2 <- tbl_data_sol12p01$PWG vec_y <- c(vec_y1, vec_y2) cat(paste(rmdhelp::bcolumn_vector(pvec = vec_y, ps_name = 'y', ps_env = '$$'), collapse = '\n'), '\n')
The right-hand side is computed as
$$ \left[ \begin{array}{c} X^TR^{-1}y \ Z^TR^{-1}y \end{array} \right] $$
The solutions are
### # coefficient matrix mat_xtrinvz <- t(mat_x) %*% mat_rinv %*% mat_z mat_ztrinvzginv <- t(mat_z) %*% mat_rinv %*% mat_z + mat_ginv mat_coeff <- rbind(cbind(mat_xtrinvx, mat_xtrinvz), cbind(t(mat_xtrinvz), mat_ztrinvzginv)) ### # right-hand side mat_rhs <- rbind(t(mat_x) %*% mat_rinv %*% vec_y, t(mat_z) %*% mat_rinv %*% vec_y) vec_sol <- solve(mat_coeff, mat_rhs) vec_beta_hat <- c("\\widehat{\\beta_{M,WWG}}", "\\widehat{\\beta_{F,WWG}}", "\\widehat{\\beta_{M,PWG}}", "\\widehat{\\beta_{F,PWG}}") vec_a1_hat <- c("\\widehat{u_{1,WWG}}", "\\widehat{u_{2,WWG}}", "\\widehat{u_{3,WWG}}", "\\widehat{u_{4,WWG}}", "\\widehat{u_{5,WWG}}", "\\widehat{u_{6,WWG}}", "\\widehat{u_{7,WWG}}", "\\widehat{u_{8,WWG}}") vec_a2_hat <- c("\\widehat{u_{1,PWG}}", "\\widehat{u_{2,PWG}}", "\\widehat{u_{3,PWG}}", "\\widehat{u_{4,PWG}}", "\\widehat{u_{5,PWG}}", "\\widehat{u_{6,PWG}}", "\\widehat{u_{7,PWG}}", "\\widehat{u_{8,PWG}}") vec_a_hat <- c(vec_a1_hat, vec_a2_hat) vec_hat_unknown <- c(vec_beta_hat, vec_a_hat) cat("$$\n") cat(paste(rmdhelp::bcolumn_vector(pvec = vec_hat_unknown), collapse = '\n'), '\n') cat(" = ") cat(paste(rmdhelp::bmatrix(pmat = round(vec_sol, digits = 4)), collapse = '\n'), '\n') cat("$$\n")
\vspace{3ex}
cat(cnt$out(ps_suffix = "Comparison of Reliabilites"), "\n")
Compare the predicted breeding values and the reliabilites obtained as results of Problem 1 with results from two univariate analyses for the same traits are used in Problem 1. All parameters can be taken from Problem 1.
Solution
For a predicted breeding value $\hat{u}_i$, the reliability $B_i$ is computed as
$$B_i = r_{u,\hat{u}}^2 = 1 - \frac{PEV(\hat{u}i)}{var(u_i)} = 1 - \frac{C{ii}^{22}}{var(u_i)}$$
where $C_{ii}^{22}$ are obtained from the inverse coefficient matrix of the mixed model equations. Just as a reminder, we can write the mixed model equations (MME) as
$$M \cdot s = r$$ with the vectors $r$ and $s$ corresponding to the right-hand side and to the unknowns of the MME. Hence
$$r = \left[\begin{array}{c}X^TR^{-1}y \ Z^TR^{-1}y \end{array}\right]$$ and
$$s = \left[\begin{array}{c}\hat{\beta} \ \hat{u} \end{array}\right]$$ The matrix $C^{22}$ is taken from the inverse coefficient matrix.
$$ M^{-1} = \left[ \begin{array}{lr} X^TR^{-1}X & X^TR^{-1}Z \ Z^TR^{-1}X & Z^TR^{-1}Z + G^{-1} \end{array} \right]^{-1} = \left[ \begin{array}{lr} C^{11} & C^{12} \ C^{21} & C^{22} \end{array} \right] $$
For the two univariate analyses, we get the solutions for the fixed effects and the breeding values and their reliabilities as follows
- WWG: estimates $s_{WWG}$ and reliabilites $B_{WWG}$
lambda1 <- mat_r0[1,1] / mat_g0[1,1] matxtx1 <- crossprod(mat_x1) matxtz1 <- crossprod(mat_x1,mat_z1) matztx1 <- t(matxtz1) matztzainv1 <- crossprod(mat_z1) + mat_ainv * lambda1 mat_coeff1 <- rbind(cbind(matxtx1, matxtz1), cbind(matztx1, matztzainv1)) mat_rhs1 <- rbind(crossprod(mat_x1, vec_y1), crossprod(mat_z1, vec_y1)) mat_sol1 <- solve(mat_coeff1, mat_rhs1)
cat(paste(rmdhelp::bmatrix(pmat = round(mat_sol1, digits = 4), ps_name = 's_{WWG}', ps_env = '$$'), collapse = '\n'), '\n')
mat_coeff_inv1 <- solve(mat_coeff1 / mat_r0[1,1]) n_nr_sol <- dim(mat_coeff_inv1)[1] vec_c22_diag1 <- diag(mat_coeff_inv1)[(n_nr_sol - n_nr_animal + 1):n_nr_sol] vec_rel1 <- 1 - vec_c22_diag1 / mat_g0[1,1]
cat(paste(rmdhelp::bcolumn_vector(pvec = round(vec_rel1, digits = 4), ps_name = 'B_{WWG}', ps_env = '$$'), collapse = '\n'), '\n')
- PWG
lambda2 <- mat_r0[2,2] / mat_g0[2,2] matxtx2 <- crossprod(mat_x2) matxtz2 <- crossprod(mat_x2,mat_z2) matztx2 <- t(matxtz2) matztzainv2 <- crossprod(mat_z2) + mat_ainv * lambda2 mat_coeff2 <- rbind(cbind(matxtx2, matxtz2), cbind(matztx2, matztzainv2)) mat_rhs2 <- rbind(crossprod(mat_x2, vec_y2), crossprod(mat_z2, vec_y2)) mat_sol2 <- solve(mat_coeff2, mat_rhs2)
cat(paste(rmdhelp::bmatrix(pmat = round(mat_sol2, digits = 4), ps_name = 's_{PWG}', ps_env = '$$'), collapse = '\n'), '\n')
mat_coeff_inv2 <- solve(mat_coeff2 / mat_r0[2,2]) n_nr_sol <- dim(mat_coeff_inv2)[1] vec_c22_diag2 <- diag(mat_coeff_inv2)[(n_nr_sol - n_nr_animal + 1):n_nr_sol] vec_rel2 <- 1 - vec_c22_diag2 / mat_g0[2,2]
cat(paste(rmdhelp::bcolumn_vector(pvec = round(vec_rel2, digits = 4), ps_name = 'B_{PWG}', ps_env = '$$'), collapse = '\n'), '\n')
The reliabilities from the bivariate analysis are obtained as
mat_coeff_inv <- solve(mat_coeff) n_nr_sol <- dim(mat_coeff_inv)[1] vec_c22_diag <- diag(mat_coeff_inv)[(n_nr_sol - n_nr_trait * n_nr_animal + 1):n_nr_sol] vec_rel <- 1 - vec_c22_diag / c(rep(mat_g0[1,1], n_nr_animal), rep(mat_g0[2,2], n_nr_animal))
cat(paste(rmdhelp::bcolumn_vector(pvec = round(vec_rel, digits = 4), ps_name = 'B', ps_env = '$$'), collapse = '\n'), '\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.