In charlotte-ngs/asmss2022: Applied Statistical Methods in Animal Science - Spring Semester 2022
knitr::opts_chunk$set(echo = TRUE)
Problem 1: Sire Model
s_ex10_p01_data_path <- "https://charlotte-ngs.github.io/asmss2022/data/asm_ped_sim_data.csv"
if (!params$isonline)
s_ex10_p01_data_path <- file.path(here::here(), "docs", "data", "asm_ped_sim_data.csv")
sigma_u2 <- 9
sigma_s2 <- sigma_u2 / 4
sigma_e2 <- 36
Use the following dataset to predict breeding values using a sire model. The dataset is available from
cat(s_ex10_p01_data_path, "\n")
Hints
- The variance component $\sigma_s^2$ of the sire effect can be assumed to be $
r sigma_s2$.
- The variance component $\sigma_e^2$ of the random resiudals is $
r sigma_e2$.
- Sex is modelled as a fixed effect.
- The sire pedigree relationship can be computed using the
pedigreemm package.
Solution
- Read the data
s_ex10_p01_data_path <- "https://charlotte-ngs.github.io/asmss2022/data/asm_ped_sim_data.csv"
tbl_ex10_p01 <- readr::read_csv(s_ex10_p01_data_path)
- Compute the inverse sire relationship matrix
vec_sire <- unique(tbl_ex10_p01$SIRE)
vec_sire <- vec_sire[!is.na(vec_sire)]
n_nr_sire <- length(vec_sire)
ped_sire <- pedigreemm::pedigree(sire = c(NA, NA, 2),
dam = rep(NA, n_nr_sire),
label = as.character(vec_sire))
mat_A_inv_sire <- as.matrix(pedigreemm::getAInv(ped = ped_sire))
mat_A_inv_sire
- Setup the mixed model equations. The mixed model equations are
$$
\left[
\begin{array}{cc}
X^TX & X^TZ \
Z^TX & Z^TZ + \lambda_s * A_s^{-1}
\end{array}
\right]
\left[
\begin{array}{c}
\hat{b} \
\hat{s}
\end{array}
\right]
=
\left[
\begin{array}{c}
X^Ty \
Z^Ty
\end{array}
\right]
$$
where $\lambda_s = \sigma_e^2 / \sigma_s^2$.
The components of the mixed model equations are shown in the following table
tbl_mme <- tibble::tibble(Component = c("$X$", "$Z$", "$y$", "$\\lambda_s$", "$A_s^{-1}$"),
Description = c("Given in the data",
"Given in the data",
"Given in the data",
"Given by variance components",
"Computed above"))
knitr::kable(tbl_mme,
booktabs = TRUE,
longtable = TRUE,
escape = FALSE)
The matrix $X$
mat_X <- model.matrix(lm(P ~ 0 + SEX, data = tbl_ex10_p01))
attr(mat_X, "assign") <- NULL
attr(mat_X, "contrasts") <- NULL
colnames(mat_X) <- NULL
mat_X
The matrix $Z$
mat_Z <- model.matrix(lm(P ~ 0 + as.factor(SIRE), data = tbl_ex10_p01))
attr(mat_Z, "assign") <- NULL
attr(mat_Z, "contrasts") <- NULL
colnames(mat_Z) <- NULL
mat_Z
The vector $y$
vec_y <- tbl_ex10_p01$P
vec_y
The mixed model equations are
mat_xtx <- crossprod(mat_X)
mat_xtz <- crossprod(mat_X, mat_Z)
mat_ztx <- t(mat_xtz)
lambda_s <- sigma_e2 / sigma_s2
mat_ztz_a_inv_lambda <- crossprod(mat_Z) + lambda_s * mat_A_inv_sire
mat_coef <- rbind(cbind(mat_xtx, mat_xtz), cbind(mat_ztx, mat_ztz_a_inv_lambda))
mat_xty <- crossprod(mat_X, vec_y)
mat_zty <- crossprod(mat_Z, vec_y)
mat_rhs <- rbind(mat_xty, mat_zty)
- Solve mixed model equations. The solution is obtained by
mat_sol_sire <- solve(mat_coef, mat_rhs)
mat_sol_sire
The solution for the fixed effects are
mat_sol_sire[1:2,]
The predicted breeding values are
mat_sol_sire[3:nrow(mat_sol_sire),]
Problem 2: Animal Model
s_ex10_p02_data_path <- "https://charlotte-ngs.github.io/asmss2022/data/asm_ped_sim_data.csv"
if (!params$isonline)
s_ex10_p02_data_path <- file.path(here::here(), "docs", "data", "asm_ped_sim_data.csv")
sigma_u2 <- 9
sigma_s2 <- sigma_u2 / 4
sigma_e2 <- 36
Use the same dataset as in Problem 1 to predict breeding values, but use an animal model instead of a sire model. The dataset is available from
cat(s_ex10_p02_data_path, "\n")
Hints
- The variance component $\sigma_u^2$ of the breeding value can be assumed to be $
r sigma_u2$.
- The variance component $\sigma_e^2$ of the random resiudals is $
r sigma_e2$.
- Sex is modelled as a fixed effect.
- The numerator relationship matrix can be computed using the
pedigreemm package.
Solution
- Read the data
s_ex10_p02_data_path <- "https://charlotte-ngs.github.io/asmss2022/data/asm_ped_sim_data.csv"
tbl_ex10_p02 <- readr::read_csv(s_ex10_p02_data_path)
- Compute the inverse sire relationship matrix
ped <- pedigreemm::pedigree(sire = c(rep(NA, 4), tbl_ex10_p02$SIRE),
dam = c(rep(NA, 4), tbl_ex10_p02$DAM),
label = as.character(c(1:4, tbl_ex10_p02$ID)))
mat_A_inv <- as.matrix(pedigreemm::getAInv(ped = ped))
mat_A_inv
- Setup the mixed model equations. The mixed model equations are
$$
\left[
\begin{array}{cc}
X^TX & X^TZ \
Z^TX & Z^TZ + \lambda * A^{-1}
\end{array}
\right]
\left[
\begin{array}{c}
\hat{b} \
\hat{u}
\end{array}
\right]
=
\left[
\begin{array}{c}
X^Ty \
Z^Ty
\end{array}
\right]
$$
where $\lambda = \sigma_e^2 / \sigma_u^2$.
The components of the mixed model equations are shown in the following table
tbl_mme <- tibble::tibble(Component = c("$X$", "$Z$", "$y$", "$\\lambda$", "$A^{-1}$"),
Description = c("Given in the data",
"Given in the data",
"Given in the data",
"Given by variance components",
"Computed above"))
knitr::kable(tbl_mme,
booktabs = TRUE,
longtable = TRUE,
escape = FALSE)
The matrix $X$
mat_X <- model.matrix(lm(P ~ 0 + SEX, data = tbl_ex10_p02))
attr(mat_X, "assign") <- NULL
attr(mat_X, "contrasts") <- NULL
colnames(mat_X) <- NULL
mat_X
The matrix $Z$
# model matrix from data
mat_Z <- model.matrix(lm(P ~ 0 + as.factor(ID), data = tbl_ex10_p02))
attr(mat_Z, "assign") <- NULL
attr(mat_Z, "contrasts") <- NULL
colnames(mat_Z) <- NULL
# add founders
mat_Z <- cbind(matrix(0, nrow = nrow(mat_Z), ncol = 4), mat_Z)
mat_Z
The vector $y$
vec_y <- tbl_ex10_p02$P
vec_y
The mixed model equations are
mat_xtx <- crossprod(mat_X)
mat_xtz <- crossprod(mat_X, mat_Z)
mat_ztx <- t(mat_xtz)
lambda <- sigma_e2 / sigma_u2
mat_ztz_a_inv_lambda <- crossprod(mat_Z) + lambda * mat_A_inv
mat_coef <- rbind(cbind(mat_xtx, mat_xtz), cbind(mat_ztx, mat_ztz_a_inv_lambda))
mat_xty <- crossprod(mat_X, vec_y)
mat_zty <- crossprod(mat_Z, vec_y)
mat_rhs <- rbind(mat_xty, mat_zty)
- Solve mixed model equations. The solution is obtained by
mat_sol <- solve(mat_coef, mat_rhs)
mat_sol
The solution for the fixed effects are
mat_sol[1:2,]
The predicted breeding values are
mat_sol[3:nrow(mat_sol),]
Problem 3: Model Comparison
Compare the order of the predicted breeding values for the sires from the sire model and from the animal model.
Solution
- Sire model
order(mat_sol_sire[3:nrow(mat_sol_sire)], decreasing = TRUE)
- Animal model
order(mat_sol[3:nrow(mat_sol)], decreasing = TRUE)
The order of the sires is the same under both models
charlotte-ngs/asmss2022 documentation built on June 7, 2022, 1:33 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 = TRUE)
Problem 1: Sire Model
s_ex10_p01_data_path <- "https://charlotte-ngs.github.io/asmss2022/data/asm_ped_sim_data.csv" if (!params$isonline) s_ex10_p01_data_path <- file.path(here::here(), "docs", "data", "asm_ped_sim_data.csv") sigma_u2 <- 9 sigma_s2 <- sigma_u2 / 4 sigma_e2 <- 36
Use the following dataset to predict breeding values using a sire model. The dataset is available from
cat(s_ex10_p01_data_path, "\n")
Hints
- The variance component $\sigma_s^2$ of the sire effect can be assumed to be $
r sigma_s2$. - The variance component $\sigma_e^2$ of the random resiudals is $
r sigma_e2$. - Sex is modelled as a fixed effect.
- The sire pedigree relationship can be computed using the
pedigreemmpackage.
Solution
- Read the data
s_ex10_p01_data_path <- "https://charlotte-ngs.github.io/asmss2022/data/asm_ped_sim_data.csv" tbl_ex10_p01 <- readr::read_csv(s_ex10_p01_data_path)
- Compute the inverse sire relationship matrix
vec_sire <- unique(tbl_ex10_p01$SIRE) vec_sire <- vec_sire[!is.na(vec_sire)] n_nr_sire <- length(vec_sire) ped_sire <- pedigreemm::pedigree(sire = c(NA, NA, 2), dam = rep(NA, n_nr_sire), label = as.character(vec_sire)) mat_A_inv_sire <- as.matrix(pedigreemm::getAInv(ped = ped_sire)) mat_A_inv_sire
- Setup the mixed model equations. The mixed model equations are
$$ \left[ \begin{array}{cc} X^TX & X^TZ \ Z^TX & Z^TZ + \lambda_s * A_s^{-1} \end{array} \right] \left[ \begin{array}{c} \hat{b} \ \hat{s} \end{array} \right] = \left[ \begin{array}{c} X^Ty \ Z^Ty \end{array} \right] $$
where $\lambda_s = \sigma_e^2 / \sigma_s^2$.
The components of the mixed model equations are shown in the following table
tbl_mme <- tibble::tibble(Component = c("$X$", "$Z$", "$y$", "$\\lambda_s$", "$A_s^{-1}$"), Description = c("Given in the data", "Given in the data", "Given in the data", "Given by variance components", "Computed above")) knitr::kable(tbl_mme, booktabs = TRUE, longtable = TRUE, escape = FALSE)
The matrix $X$
mat_X <- model.matrix(lm(P ~ 0 + SEX, data = tbl_ex10_p01)) attr(mat_X, "assign") <- NULL attr(mat_X, "contrasts") <- NULL colnames(mat_X) <- NULL mat_X
The matrix $Z$
mat_Z <- model.matrix(lm(P ~ 0 + as.factor(SIRE), data = tbl_ex10_p01)) attr(mat_Z, "assign") <- NULL attr(mat_Z, "contrasts") <- NULL colnames(mat_Z) <- NULL mat_Z
The vector $y$
vec_y <- tbl_ex10_p01$P vec_y
The mixed model equations are
mat_xtx <- crossprod(mat_X) mat_xtz <- crossprod(mat_X, mat_Z) mat_ztx <- t(mat_xtz) lambda_s <- sigma_e2 / sigma_s2 mat_ztz_a_inv_lambda <- crossprod(mat_Z) + lambda_s * mat_A_inv_sire mat_coef <- rbind(cbind(mat_xtx, mat_xtz), cbind(mat_ztx, mat_ztz_a_inv_lambda)) mat_xty <- crossprod(mat_X, vec_y) mat_zty <- crossprod(mat_Z, vec_y) mat_rhs <- rbind(mat_xty, mat_zty)
- Solve mixed model equations. The solution is obtained by
mat_sol_sire <- solve(mat_coef, mat_rhs) mat_sol_sire
The solution for the fixed effects are
mat_sol_sire[1:2,]
The predicted breeding values are
mat_sol_sire[3:nrow(mat_sol_sire),]
Problem 2: Animal Model
s_ex10_p02_data_path <- "https://charlotte-ngs.github.io/asmss2022/data/asm_ped_sim_data.csv" if (!params$isonline) s_ex10_p02_data_path <- file.path(here::here(), "docs", "data", "asm_ped_sim_data.csv") sigma_u2 <- 9 sigma_s2 <- sigma_u2 / 4 sigma_e2 <- 36
Use the same dataset as in Problem 1 to predict breeding values, but use an animal model instead of a sire model. The dataset is available from
cat(s_ex10_p02_data_path, "\n")
Hints
- The variance component $\sigma_u^2$ of the breeding value can be assumed to be $
r sigma_u2$. - The variance component $\sigma_e^2$ of the random resiudals is $
r sigma_e2$. - Sex is modelled as a fixed effect.
- The numerator relationship matrix can be computed using the
pedigreemmpackage.
Solution
- Read the data
s_ex10_p02_data_path <- "https://charlotte-ngs.github.io/asmss2022/data/asm_ped_sim_data.csv" tbl_ex10_p02 <- readr::read_csv(s_ex10_p02_data_path)
- Compute the inverse sire relationship matrix
ped <- pedigreemm::pedigree(sire = c(rep(NA, 4), tbl_ex10_p02$SIRE), dam = c(rep(NA, 4), tbl_ex10_p02$DAM), label = as.character(c(1:4, tbl_ex10_p02$ID))) mat_A_inv <- as.matrix(pedigreemm::getAInv(ped = ped)) mat_A_inv
- Setup the mixed model equations. The mixed model equations are
$$ \left[ \begin{array}{cc} X^TX & X^TZ \ Z^TX & Z^TZ + \lambda * A^{-1} \end{array} \right] \left[ \begin{array}{c} \hat{b} \ \hat{u} \end{array} \right] = \left[ \begin{array}{c} X^Ty \ Z^Ty \end{array} \right] $$
where $\lambda = \sigma_e^2 / \sigma_u^2$.
The components of the mixed model equations are shown in the following table
tbl_mme <- tibble::tibble(Component = c("$X$", "$Z$", "$y$", "$\\lambda$", "$A^{-1}$"), Description = c("Given in the data", "Given in the data", "Given in the data", "Given by variance components", "Computed above")) knitr::kable(tbl_mme, booktabs = TRUE, longtable = TRUE, escape = FALSE)
The matrix $X$
mat_X <- model.matrix(lm(P ~ 0 + SEX, data = tbl_ex10_p02)) attr(mat_X, "assign") <- NULL attr(mat_X, "contrasts") <- NULL colnames(mat_X) <- NULL mat_X
The matrix $Z$
# model matrix from data mat_Z <- model.matrix(lm(P ~ 0 + as.factor(ID), data = tbl_ex10_p02)) attr(mat_Z, "assign") <- NULL attr(mat_Z, "contrasts") <- NULL colnames(mat_Z) <- NULL # add founders mat_Z <- cbind(matrix(0, nrow = nrow(mat_Z), ncol = 4), mat_Z) mat_Z
The vector $y$
vec_y <- tbl_ex10_p02$P vec_y
The mixed model equations are
mat_xtx <- crossprod(mat_X) mat_xtz <- crossprod(mat_X, mat_Z) mat_ztx <- t(mat_xtz) lambda <- sigma_e2 / sigma_u2 mat_ztz_a_inv_lambda <- crossprod(mat_Z) + lambda * mat_A_inv mat_coef <- rbind(cbind(mat_xtx, mat_xtz), cbind(mat_ztx, mat_ztz_a_inv_lambda)) mat_xty <- crossprod(mat_X, vec_y) mat_zty <- crossprod(mat_Z, vec_y) mat_rhs <- rbind(mat_xty, mat_zty)
- Solve mixed model equations. The solution is obtained by
mat_sol <- solve(mat_coef, mat_rhs) mat_sol
The solution for the fixed effects are
mat_sol[1:2,]
The predicted breeding values are
mat_sol[3:nrow(mat_sol),]
Problem 3: Model Comparison
Compare the order of the predicted breeding values for the sires from the sire model and from the animal model.
Solution
- Sire model
order(mat_sol_sire[3:nrow(mat_sol_sire)], decreasing = TRUE)
- Animal model
order(mat_sol[3:nrow(mat_sol)], decreasing = TRUE)
The order of the sires is the same under both models
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.