In charlotte-ngs/lbgfs2021: Livestock Breeding and Genomics (ETHZ) - Fall Semester 2021
knitr::opts_chunk$set(echo = TRUE)
knitr::knit_hooks$set(hook_convert_odg = rmdhelp::hook_convert_odg)
library(dplyr)
cnt <- rmdhelp::R6ClassCount$new()
cnt$set_prefix(ps_prefix = "## Problem")
# Assign Points for Q1
lPointsQ1 <- list(TaskA = 20,
TaskB = 20,
TaskC = 5,
TaskD = 0)
nPointQ1Total <- sum(unlist(lPointsQ1))
# Assign Points for Q2
lPointsQ2 <- list(TaskA = 10,
TaskB = 30,
TaskC = 0)
nPointQ2Total <- sum(unlist(lPointsQ2))
# Assign Points for Q3
lPointsQ3 <- list(TaskA = 12,
TaskB = 6,
TaskC = 0)
nPointQ3Total <- sum(unlist(lPointsQ3))
# Assign Points for Q4
lPointsQ4 <- list(TaskA = 4,
TaskB = 10,
TaskC = 2)
nPointQ4Total <- sum(unlist(lPointsQ4))
# Assign Points for Q5
lPointsQ5 <- list(TaskA = 5,
TaskB = 5,
TaskC = 0)
nPointQ5Total <- sum(unlist(lPointsQ5))
# compute overal sum of points
nPointOverallTotal <- nPointQ1Total + nPointQ2Total + nPointQ3Total + nPointQ4Total + nPointQ5Total
\thispagestyle{empty}
\begin{tabular}{l}
ETH Zurich \
D-USYS\
Institute of Agricultural Sciences\
\end{tabular}
\vspace{15ex}
\begin{center}
\huge
Solutions To Exam\ \vspace{1ex}
Livestock Breeding and Genomics \ \vspace{1ex}
FS 2020 \
\normalsize
\vspace{7ex}
Peter von Rohr
\end{center}
\vspace{7ex}
\begin{tabular}{p{5cm}lr}
& \textsc{Date} & \textsc{\emph{18. December 2020}} \
& \textsc{Begin} & \textsc{\emph{09:15 }}\
& \textsc{End} & \textsc{\emph{11:15 }}\
\end{tabular}
\vspace{5ex}
\large
\begin{tabular}{p{2.5cm}p{3cm}p{6cm}}
& Name: & \
& & \
& Legi-Nr: & \
\end{tabular}
\normalsize
\vspace{9ex}
\begin{center}
\begin{tabular}{|p{3cm}|c|c|}
\hline
Problem & Maximum Number of Points & Number of Points Reached \
\hline
1 & r nPointQ1Total & \
\hline
2 & r nPointQ2Total & \
\hline
3 & r nPointQ3Total & \
\hline
4 & r nPointQ4Total & \
\hline
5 & r nPointQ5Total & \
\hline
Total & r nPointOverallTotal & \
\hline
\end{tabular}
\end{center}
\clearpage
\pagebreak
tbl_ped_nrm <- tibble::tibble(ID = c(5,6,7,8,9,10), Sex = c('M','F','M','F','M','M'), Sire = c(3,1,5,5,7,9), Dam = c(2,4,6,6,8,8))
dam_id <- tbl_ped_nrm$Dam[length(tbl_ped_nrm$Dam)]
sire_id <- tbl_ped_nrm$Sire[length(tbl_ped_nrm$Sire)]
cat(cnt$out(ps_suffix = "Numerator Relationship Matrix and Inbreeding"), "\n")
We are given the following pedigree.
\textit{Das folgende Pedigree ist gegeben.}
\begin{center}
knitr::kable(tbl_ped_nrm,
booktabs = TRUE,
escape = FALSE,
format = 'latex')
\end{center}
\begin{enumerate}
\item[a)] Compute the additive genetic relationship matrix $A$ for the above pedigree.
\textit{Berechnen Sie die additiv-genetische Verwandtschaftsmatrix $A$ für das oben angegebene Pedigree}
\points{r lPointsQ1$TaskA}
\end{enumerate}
\solstart
The numerator relationship matrix is computed using pedigreemm::getA(). In a first step, we have to extend the pedigree to contain the founder animals in the ID-column.
vec_founder_sire <- setdiff(tbl_ped_nrm$Sire, tbl_ped_nrm$ID)
n_nr_founder_sire <- length(vec_founder_sire)
vec_founder_dam <- setdiff(tbl_ped_nrm$Dam, tbl_ped_nrm$ID)
n_nr_founder_dam <- length(vec_founder_dam)
# check that founder_sire and founder_dam are not the same
if (length(intersect(vec_founder_sire, vec_founder_dam)) != 0)
stop(" * ERROR: Founder sires and founder dams are not exclusive")
tbl_ped_nrm_ext <- dplyr::bind_rows(
tibble::tibble(ID = vec_founder_sire[order(vec_founder_sire)],
Sex = rep('M', n_nr_founder_sire),
Sire = rep(NA, n_nr_founder_sire),
Dam = rep(NA, n_nr_founder_sire)),
tibble::tibble(ID = vec_founder_dam[order(vec_founder_dam)],
Sex = rep('F', n_nr_founder_dam),
Sire = rep(NA, n_nr_founder_dam),
Dam = rep(NA, n_nr_founder_dam)),
tbl_ped_nrm)
ped_nrm <- pedigreemm::pedigree(sire = tbl_ped_nrm_ext$Sire,
dam = tbl_ped_nrm_ext$Dam,
label = as.character(tbl_ped_nrm_ext$ID))
(matA_nrm <- as.matrix(pedigreemm::getA(ped = ped_nrm)))
\solend
\clearpage
\pagebreak
\begin{enumerate}
\item[b)] Compute the inbreeding coefficients of all animals in the given pedigree. Complete the following table and indicate which of the animals are inbred.
\textit{Berechnen Sie die Inzuchtkoeffizienten aller Tiere im gegebenen Pedigree. Vervollständigen Sie die folgende Tabelle und geben an, welche Tiere ingezüchtet sind.}
\points{r lPointsQ1$TaskB}
\end{enumerate}
\vspace{3ex}
n_nr_animal <- max(tbl_ped_nrm$ID)
tbl_inb_empty <- tibble::tibble(ID = 1:n_nr_animal,
`Inbreeding Coefficient` = rep('', n_nr_animal),
`Animal Inbred (TRUE/FALSE)` = rep('', n_nr_animal))
knitr::kable(tbl_inb_empty,
booktabs = TRUE,
escape = FALSE,
format = 'latex')
\solstart
\vspace{5ex}
The numeric solution is
vec_inbr_yn <- ifelse(diag(matA_nrm) > 1, "TRUE", "FALSE")
vec_inbr_coef <- diag(matA_nrm) - 1
# tibble for table
(tbl_inb_sol <- tibble::tibble(ID = tbl_ped_nrm_ext$ID,
`Inbreeding Coefficient` = round(vec_inbr_coef, digits = 3),
`Animal Inbred (TRUE/FALSE)` = vec_inbr_yn))
\solend
\clearpage
\pagebreak
\begin{enumerate}
\item[c)] Assume that dam r dam_id and sire r sire_id are mated. What is the inbreeding coefficient of their offspring?
\textit{Wir nehmen an, dass die Mutter r dam_id mit dem Vater r sire_id angepaart wird. Wie gross ist der Inzuchtkoeffizient des Nachkommens aus dieser Paarung?}
\points{r lPointsQ1$TaskC}
\end{enumerate}
\solstart
The inbreeding coefficient of the offspring of dam r dam_id and sire r sire_id is half of the relationship between the parents.
0.5 * matA_nrm[sire_id, dam_id]
\solend
\clearpage
\pagebreak
tbl_data_pbv <- tibble::tibble(ID = c(4,5,6,7,8),
Sire = c(1,3,4,4,6),
Dam = c(2,2,5,5,7),
Herd = c('Planta','Moos','Moos','Moos','Moos'),
`Phenotypic Observation` = c(7.53,8.48,0.26,6.6,2.44))
n_var_u <- 7.2
n_var_e <- 28.8
n_h2 <- 0.2
cat(cnt$out(ps_suffix = "Prediction of Breeding Values"), "\n")
The following dataset is used to predict breeding values.
\textit{Der folgende Datensatz soll für die Schätzung von Zuchtwerten verwendet werden.}
knitr::kable(tbl_data_pbv,
booktabs = TRUE,
escape = FALSE,
format = 'latex')
Assumptions
- the error variance is: $\sigma_e^2 =
r n_var_e$
- the heritability is: $h^2 =
r n_h2$
- the genetic additive variance is: $\sigma_u^2 =
r n_var_u$
- the population mean $\mu$ can be taken as the mean of the given phenotypic observations
\textit{Annahmen}
- \textit{die Restvarianz beträgt:} $\sigma_e^2 =
r n_var_e$
- \textit{die Heritabilität beträgt:} $h^2 =
r n_h2$
- \textit{die genetisch-additive Varianz beträgt:} $\sigma_u^2 =
r n_var_u$
- \textit{das Populationsmittel $\mu$ kann als Mittelwert der phänotypischen Beobachtungen angenommen werden.}
\vspace{3ex}
\begin{enumerate}
\item[a)] Use the regression method to predict breeding values based on own-performance records for the animals in the table given above.
\textit{Verwenden Sie die Regressionsmethode zur Schätzung der Zuchtwerte basierend auf der Eigenleistung der Tiere in der obigen Tabelle.}
\points{r lPointsQ2$TaskA}
\end{enumerate}
\clearpage
\pagebreak
\solstart
According to the regression method, the predicted breeding value ($\hat{u}_i$) for animal $i$ is
$$\hat{u}_i = h^2(y_i - \mu)$$
where $y_i$ is the phenotypic observation of animal $i$, $h^2$ is the heritability and $\mu$ is the population mean.
\vspace{3ex}
vec_y <- tbl_data_pbv$`Phenotypic Observation`
n_mu <- mean(vec_y)
vec_uhat <- n_h2 * (vec_y - n_mu)
tbl_sol <- tibble::tibble(ID = tbl_data_pbv$ID, PBV = vec_uhat)
knitr::kable(tbl_sol,
booktabs = TRUE,
escape = FALSE,
format = 'latex')
\solend
\clearpage
\pagebreak
\begin{enumerate}
\item[b)] Use a BLUP animal model to predict the breeding values for all animals in the pedigree based on the data given in the table above. Specify the model to predict breeding values, name all model components, compute the expected values and the variance-covariance matrices for all random model components. Insert the information from the above table into the model components where possible. Set up the mixed model equations and compute the solutions for the estimates of fixed effects and for the predicted breeding values.
\textit{Verwenden Sie ein BLUP Tiermodell zur Schätzung der Zuchtwerte aller Tiere im Pedigree basierend auf den Daten in der obigen Tabelle. Spezifizieren Sie das Modell für die Schätzung der Zuchtwerte, benennen Sie alle Modellkomponenten, berechnen Sie die Erwartungswerte und die Varianz-Kovarianz Matrizen aller zufälligen Effekte im Modell. Setzen Sie die verfügbaren Information aus dem Datensatz in die Modellkomponenten ein. Stellen Sie die Mischmodellgleichungen auf und berechnen Sie die Schätzungen der fixen Effekte und der Zuchtwerte.}
\points{r lPointsQ2$TaskB}
\end{enumerate}
\solstart
Model
The BLUP animal model corresponds to the following mixed-effects model
$$y = Xb + Zu + e$$
where $y$ is the vector of observations, $b$ is the vector of fixed effects for the two herds, $u$ is the vector of breeding values for all animals in the pedigree and $e$ is the vector of random residuals. The design matrices $X$ and $Z$ link the fixed effects and the breeding values to the observations, respectively.
Expected Values and Variance-Covariance Matrices
Expected values and variance-covariance matrices of the random components $y$, $u$ and $e$
$$
E \left[\begin{array}{c} y \ u\ e\end{array}\right] = \left[\begin{array}{c} Xb \ 0 \ 0\end{array}\right] \text{, }
var \left[\begin{array}{c} y \ u \ e\end{array}\right] = \left[\begin{array}{ccc} ZGZ^T+R & ZG & 0 \ GZ^T & G & 0 \ 0 & 0 & R \end{array}\right]
$$
Informations in Model Components
The model vectors are
vec_y <- tbl_data_pbv$`Phenotypic Observation`
n_nr_obs <- length(vec_y)
n_nr_fix <- nlevels(as.factor(tbl_data_pbv$Herd))
vec_b <- rmdhelp::vecGetVecElem(psBaseElement = 'b', pnVecLen = n_nr_fix, psResult = 'latex')
n_nr_animal_pedigree <- length(unique(c(tbl_data_pbv$ID, tbl_data_pbv$Sire, tbl_data_pbv$Dam)))
vec_u <- rmdhelp::vecGetVecElem(psBaseElement = 'u', pnVecLen = n_nr_animal_pedigree, psResult = 'latex')
vec_e <- rmdhelp::vecGetVecElem(psBaseElement = 'e', pnVecLen = n_nr_obs, psResult = 'latex')
cat('$$\n')
cat(paste(rmdhelp::bcolumn_vector(pvec = vec_y, ps_name = 'y'), collapse = '\n'))
cat('\\text{, }')
cat(paste(rmdhelp::bcolumn_vector(pvec = vec_b, ps_name = 'b'), collapse = '\n'))
cat('\\text{, }')
cat(paste(rmdhelp::bcolumn_vector(pvec = vec_u, ps_name = 'u'), collapse = '\n'))
cat('\\text{, }')
cat(paste(rmdhelp::bcolumn_vector(pvec = vec_e, ps_name = 'e'), collapse = '\n'), '\n')
cat('$$\n')
The model matrices
# matrix X
model_mat <- model.matrix(lm(`Phenotypic Observation` ~ Herd, data = tbl_data_pbv))
mat_X <- model_mat[,]
dimnames(mat_X) <- NULL
mat_X[, 1] <- mat_X[, 1] - mat_X[, 2]
# matrix Z
mat_Z <- cbind(matrix(0, nrow = n_nr_obs, ncol = (n_nr_animal_pedigree-n_nr_obs)), diag(1, nrow = n_nr_obs, ncol = n_nr_obs))
# output
cat('$$\n')
cat(paste(rmdhelp::bmatrix(pmat = mat_X, ps_name = 'X'), collapse = '\n'), '\n')
cat('\\text{, }')
cat(paste(rmdhelp::bmatrix(pmat = mat_Z, ps_name = 'Z'), collapse = '\n'), '\n')
cat('$$\n')
Mixed Model Equations
\begin{equation}
\left[
\begin{array}{lr}
X^TX & X^TZ \
Z^TX & Z^TZ + A^{-1}* \lambda
\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]
\notag
\end{equation}
Solution
# nrm A
tbl_ped_pbv_ext <- rmdexam::extend_pedigree(ptbl_ped =
tbl_data_pbv[,c('ID', 'Sire', 'Dam')])
ped_pbv <- pedigreemm::pedigree(sire = tbl_ped_pbv_ext$Sire,
dam = tbl_ped_pbv_ext$Dam,
label = as.character(tbl_ped_pbv_ext$ID))
mat_Ainv_pbv <- as.matrix(pedigreemm::getAInv(ped = ped_pbv))
lambda <- n_var_e / n_var_u
# coefficient matrix
mat_xtx <- crossprod(mat_X)
mat_xtz <- crossprod(mat_X, mat_Z)
mat_ztx <- crossprod(mat_Z, mat_X)
mat_ztz_ainv_lambda <- crossprod(mat_Z) + lambda * mat_Ainv_pbv
mat_coef <- rbind(cbind(mat_xtx,mat_xtz), cbind(mat_ztx,mat_ztz_ainv_lambda))
# rhs
mat_xty <- crossprod(mat_X, vec_y)
mat_zty <- crossprod(mat_Z, vec_y)
mat_rhs <- rbind(mat_xty, mat_zty)
# solutions
(mat_sol <- solve(mat_coef, mat_rhs))
\solend
\clearpage
\pagebreak
tbl_data_gnm <- tibble::tibble(Animal = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15), `SNPLGH` = c(-1,1,1,0,-1,0,-1,1,1,0,1,0,1,-1,1), `SNPFS2` = c(0,1,1,0,0,-1,-1,0,0,-1,-1,0,0,1,0), Observation = c(12.7,46,32.7,19.3,14.8,6.7,2.4,33,29.4,4.9,14.4,19.5,25.6,19.1,25.6))
n_nr_snp <- 2
cat(cnt$out(ps_suffix = "Genomics"), "\n")
Given is the following data set of SNP-Genotyping results.
\textit{Gegeben sind die Genotypisierungsresultate in der nachfolgenden Tabelle.}
knitr::kable(tbl_data_gnm,
booktabs = TRUE,
escape = FALSE,
format = 'latex')
\vspace{3ex}
\begin{enumerate}
\item[a)] Use a marker effect model to estimate the fixed effects of both markers on the observation. Please specify the fixed-effect model that you use, name all the model components and insert the information from the data into the components where possible.
\textit{Verwenden Sie eine Marker-Effekt Modell zur Schätzung der fixen Effekte der beiden Marker auf die Beobachtung. Bitte spezifizieren Sie das fixe Modell, benennen Sie alle Modellkomponenten und fügen Sie die Daten aus der obigen Tabelle in die Modellkomponenten ein, wo das möglich ist.}
\points{r lPointsQ3$TaskA}
\end{enumerate}
\solstart
Model
The fixed effect model to estimate the marker effects is given by
$$y = Xb + e$$
where $y$ is the vector of observations, $b$ is the vector of fixed effects containing intercept and effect of both markers, $e$ is the vector of random error terms.
Information in Model
The model vectors are
vec_y <- tbl_data_gnm$Observation
n_nr_obs <- length(vec_y)
vec_b <- rmdhelp::vecGetVecElem(psBaseElement = 'b', pnVecLen = n_nr_snp, psResult = 'latex')
vec_e <- rmdhelp::vecGetVecElem(psBaseElement = 'e', pnVecLen = n_nr_obs, psResult = 'latex')
cat("$$\n")
cat(paste(rmdhelp::bcolumn_vector(pvec = vec_y, ps_name = 'y'), collapse = '\n'), '\n')
cat('\\text{, }')
cat(paste(rmdhelp::bcolumn_vector(pvec = vec_b, ps_name = 'b'), collapse = '\n'), '\n')
cat('\\text{, }')
cat(paste(rmdhelp::bcolumn_vector(pvec = vec_e, ps_name = 'e'), collapse = '\n'), '\n')
cat("$$\n")
Marker Effects
The marker effects are estimated using lm()
vec_mrk_names <- grep('SNP', colnames(tbl_data_gnm), value = TRUE)
lm_mrk_eff <- lm(as.formula(paste0('Observation ~ ',
paste0(vec_mrk_names, collapse = ' + '), collapse = '')),
data = tbl_data_gnm)
summary(lm_mrk_eff)
\solend
\clearpage
\pagebreak
\begin{enumerate}
\item[b)] Predict the direct genomic breeding values for all animals of the dataset using the marker effects estimated in Task a).
\textit{ Schätzen Sie die direkt-genomischen Zuchtwerte für alle Tiere im Datensatz unter Verwendung der aus Aufgabe a) geschätzten Markereffekte.}
\points{r lPointsQ3$TaskB}
\end{enumerate}
\solstart
Direct genomic breeding values $\hat{g}_i$ for animal $i$ corresponds to the sums of the marker effects, hence
$$\hat{g}_i = x_i^T \cdot \hat{b}$$
where $\hat{b}$ is the vector of estimated marker effects and $x_i^T$ is the row $i$ of the design matrix $X$ corresponding to animal $i$.
vec_mrk_eff <- coefficients(lm_mrk_eff)
n_intercept_pos <- which(names(vec_mrk_eff) == '(Intercept)')
if (length(n_intercept_pos) > 0) vec_mrk_eff <- vec_mrk_eff[-n_intercept_pos]
vec_snp_names <- names(vec_mrk_eff)
mat_geno_gnm <- as.matrix(tbl_data_gnm[, vec_snp_names])
(vec_hat_g <- crossprod(t(mat_geno_gnm), vec_mrk_eff))
\solend
\clearpage
\pagebreak
n_nr_base_animals <- 280
n_nr_gen <- 16
n_hom_diff <- 50
n_maf <- 0.2
n_level_inb <- 0.01
n_inb_dep <- 0.2
n_dom_dev <- 30
cat(cnt$out(ps_suffix = "Variance and Inbreeding"), "\n")
In dairy cattle semen and embryos of the breeds Brown Swiss and Holstein are often imported from North America. For the Brown Swiss breed, the North American population is based on r n_nr_base_animals female animals. The following assumptions can be made.
- Although, cattle does not follow the same mode of inheritance as the organism shown in the lecture, the computations as shown in the lecture can be used as an approximation.
- In contrast to the size $N$ of the subpopulations that was assumed to be the number of individuals, here $N$ is the number of female animals in a given subpopulation.
\textit{Samen und Embryos der Rassen Brown Swiss und Holstein werden oft aus Nordamerika importiert. Für die Rasse Brown Swiss basiert die Nordamerikanische Population auf r n_nr_base_animals weiblichen Tieren. Die folgenden Annahmen können getroffen werden.}
- \textit{Obwohl das Rind nicht den gleichen Vererbungsmodus zeigt, wie die Organismen, welche in den Vorlesungsunterlagen verwendet wurden, können die dort eingeführten Berechnungen der Inzucht als Annäherungen verwendet werden.}
- \textit{Die Populationsgrösse $N$ hier entspricht der Anzahl weiblichen Tiere.}
\vspace{3ex}
\begin{enumerate}
\item[a)] What is the expected level of inbreeding ($F$) of imported semen in the Brown Swiss breed after r n_nr_gen generations?
\textit{Welches ist der erwartet Wert der Inzucht ($F$) von importiertem Samen in der Rasse Brown Swiss nach r n_nr_gen Generationen?}
\points{r lPointsQ4$TaskA}
\end{enumerate}
\solstart
The level of inbreeding $F_t$ after $t$ generations is computed as
$$F_t = 1 - (1 - \Delta F)^t$$
with $\Delta F = 1/(2N)$ and $N$ being the number of female animals in the base population.
(inb_level_ft <- 1 - (1 - 1/(2*n_nr_base_animals))^n_nr_gen)
\solend
\clearpage
\pagebreak
\begin{enumerate}
\item[b)] Compute the between-line, the within-line and the total genetic variance for a single additive Locus where the difference between the homozygous genotypes is r n_hom_diff, the allele frequency $p = r n_maf$ and the level of inbreeding is r n_level_inb.
\textit{Berechnen Sie die Innerhalb-Linie, Zwischen-Linie und die totale genetische Varianz für einen additiven Lokus, wobei die Differenz der homozygoten Genotypen r n_hom_diff entspricht, die Allelefrequenz $p = r n_maf$ ist und der Inzuchtwert r n_level_inb ist.}
\points{r lPointsQ4$TaskB}
\end{enumerate}
\solstart
The variance according to the sources are computed as
tbl_gen_anova <- tibble::tibble(Source = c("Between lines", "Within lines","Total"),
Variance = c("$2FV_G$", "$(1-F)V_G$", "$(1+F)V_G$"))
knitr::kable(tbl_gen_anova,
booktabs = TRUE,
escape = FALSE,
format = 'latex')
n_vg <- 2 * n_maf * (1-n_maf) * (n_hom_diff/2)^2
with
$$V_G = 2pqa^2 = 2 * r n_maf * r 1-n_maf * r (n_hom_diff/2)^2 = r n_vg$$
Inserting the values, we get
tbl_gen_anova_result <- dplyr::bind_cols(tbl_gen_anova, tibble::tibble(Results = c(2*n_level_inb*n_vg,
(1-n_level_inb)*n_vg,
(1+n_level_inb)*n_vg)))
knitr::kable(tbl_gen_anova_result,
booktabs = TRUE,
escape = FALSE,
format = 'latex')
\solend
\clearpage
\pagebreak
\begin{enumerate}
\item[c)] After how many generations is the expected inbreeding depression bigger than r n_inb_dep in a population with $N = r n_nr_base_animals$ animals. The following assumptions can be made
\textit{Nach wie vielen Generationen ist die erwartete Inzuchtdepression grösser als r n_inb_dep in einer Population von $N = r n_nr_base_animals$ Tieren. Die folgenden Annahmen können getroffen werden.}
\points{r lPointsQ4$TaskC}
\end{enumerate}
- single bi-allelic locus
- minor allele frequency $p =
r n_maf$
-
dominance deviation $d = r n_dom_dev$
-
\textit{einzelner Locus mit zwei Allelen}
- \textit{Frequenz des seltenen Alleles $p =
r n_maf$}
- \textit{Dominanzabweichung $d =
r n_dom_dev$}
\solstart
n_inb_f <- n_inb_dep / (2*n_dom_dev*n_maf*(1-n_maf))
n_nr_dgen <- log(1-n_inb_f) / log(1 - 1/(2*n_nr_base_animals))
Inbreeding depression is computed as
$$M_0 - M_F = 2dp(1-p)F = r n_inb_dep$$
hence
$$F = \frac{M_0 - M_F}{2dp(1-p)} = \frac{r n_inb_dep}{2 * r n_dom_dev * r n_maf * (1- r n_maf)} = r n_inb_f$$
From
$$F_t = 1 - (1 - \Delta F)^t$$
we get
$$(1 - \Delta F)^t = 1-F_t$$
and
$$t = \frac{log(1-F_t)}{log(1 - \Delta F)} = r n_nr_dgen$$
Hence after r ceiling(n_nr_dgen) generations the inbreeding depression is bigger than r n_inb_dep
\solend
\clearpage
\pagebreak
tbl_data_anova <- tibble::tibble(Animal = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20), Farm = c(3,1,2,3,1,2,3,1,2,1,2,1,2,2,1,1,1,2,3,3), LiveWeight = c(613,621,630,614,629,611,612,614,606,621,621,623,608,603,589,599,610,595,612,616))
cat(cnt$out(ps_suffix = "Variance Components"), "\n")
We are given the following dataset for the trait live weight (LiveWeight) for cattle.
\textit{Der folgende Datensatz umfasst das Merkmal Lebendgewicht (LiveWeight) von Rindern.}
knitr::kable(tbl_data_anova,
booktabs = TRUE,
escape = FALSE,
format = 'latex')
\begin{enumerate}
\item[a)] Compute the estimate of the error variance $\sigma_e^2$ from the residuals of the fixed linear model specified below.
\textit{Schätzen Sie die Fehlervarianz $\sigma_e^2$ basierend auf den Residuen des folgenden fixen Modells.}
\points{r lPointsQ5$TaskA}
\end{enumerate}
The fixed linear model that is used is
$$y = Xb + e$$
where $y$ is the vector of all live weight values, $b$ is the vector of the effects caused by the different farms. The fixed linear model is specified in R by
tbl_data_anova$Farm <- as.factor(tbl_data_anova$Farm)
lm_lweight <- lm(LiveWeight ~ 0 + Farm, data = tbl_data_anova)
The resulting effects of the farms are
(vec_coef_lweight <- coefficients(lm_lweight))
\solstart
The esimate of the error variance is computed based on the resiudals. The residuals can be obtained by the function residuals() in R.
vec_res <- residuals(lm_lweight)
ssq_res <- crossprod(vec_res)
(n_est_res_var <- ssq_res / (nrow(tbl_data_anova)-length(vec_coef_lweight)))
The error standard deviation is
(n_est_res_sd <- sqrt(n_est_res_var))
\solend
\clearpage
\pagebreak
\begin{enumerate}
\item[b)] Verify your result from task a) with the output of the summary()-function applied to the result of the lm()-function
\textit{Verifizieren Sie das Resultat aus Aufgabe a) anhand des Outputs der summary-Funktion angewendet auf das Resultat der lm()-Funktion}
\points{r lPointsQ5$TaskB}
\end{enumerate}
\solstart
From the task, we have the result object of the lm()-function which is called lm_lweight. Applying the summary()-method leads to
summary(lm_lweight)
The number next to Residual standard error: corresponds to the estimated value of the error standard deviation.
\solend
charlotte-ngs/lbgfs2021 documentation built on Dec. 19, 2021, 3:01 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) knitr::knit_hooks$set(hook_convert_odg = rmdhelp::hook_convert_odg) library(dplyr)
cnt <- rmdhelp::R6ClassCount$new() cnt$set_prefix(ps_prefix = "## Problem")
# Assign Points for Q1 lPointsQ1 <- list(TaskA = 20, TaskB = 20, TaskC = 5, TaskD = 0) nPointQ1Total <- sum(unlist(lPointsQ1)) # Assign Points for Q2 lPointsQ2 <- list(TaskA = 10, TaskB = 30, TaskC = 0) nPointQ2Total <- sum(unlist(lPointsQ2)) # Assign Points for Q3 lPointsQ3 <- list(TaskA = 12, TaskB = 6, TaskC = 0) nPointQ3Total <- sum(unlist(lPointsQ3)) # Assign Points for Q4 lPointsQ4 <- list(TaskA = 4, TaskB = 10, TaskC = 2) nPointQ4Total <- sum(unlist(lPointsQ4)) # Assign Points for Q5 lPointsQ5 <- list(TaskA = 5, TaskB = 5, TaskC = 0) nPointQ5Total <- sum(unlist(lPointsQ5)) # compute overal sum of points nPointOverallTotal <- nPointQ1Total + nPointQ2Total + nPointQ3Total + nPointQ4Total + nPointQ5Total
\thispagestyle{empty}
\begin{tabular}{l} ETH Zurich \ D-USYS\ Institute of Agricultural Sciences\ \end{tabular}
\vspace{15ex} \begin{center} \huge Solutions To Exam\ \vspace{1ex} Livestock Breeding and Genomics \ \vspace{1ex} FS 2020 \
\normalsize \vspace{7ex} Peter von Rohr \end{center}
\vspace{7ex} \begin{tabular}{p{5cm}lr} & \textsc{Date} & \textsc{\emph{18. December 2020}} \ & \textsc{Begin} & \textsc{\emph{09:15 }}\ & \textsc{End} & \textsc{\emph{11:15 }}\ \end{tabular}
\vspace{5ex} \large \begin{tabular}{p{2.5cm}p{3cm}p{6cm}} & Name: & \ & & \ & Legi-Nr: & \ \end{tabular} \normalsize
\vspace{9ex}
\begin{center}
\begin{tabular}{|p{3cm}|c|c|}
\hline
Problem & Maximum Number of Points & Number of Points Reached \
\hline
1 & r nPointQ1Total & \
\hline
2 & r nPointQ2Total & \
\hline
3 & r nPointQ3Total & \
\hline
4 & r nPointQ4Total & \
\hline
5 & r nPointQ5Total & \
\hline
Total & r nPointOverallTotal & \
\hline
\end{tabular}
\end{center}
\clearpage \pagebreak
tbl_ped_nrm <- tibble::tibble(ID = c(5,6,7,8,9,10), Sex = c('M','F','M','F','M','M'), Sire = c(3,1,5,5,7,9), Dam = c(2,4,6,6,8,8)) dam_id <- tbl_ped_nrm$Dam[length(tbl_ped_nrm$Dam)] sire_id <- tbl_ped_nrm$Sire[length(tbl_ped_nrm$Sire)]
cat(cnt$out(ps_suffix = "Numerator Relationship Matrix and Inbreeding"), "\n")
We are given the following pedigree.
\textit{Das folgende Pedigree ist gegeben.}
\begin{center}
knitr::kable(tbl_ped_nrm, booktabs = TRUE, escape = FALSE, format = 'latex')
\end{center}
\begin{enumerate} \item[a)] Compute the additive genetic relationship matrix $A$ for the above pedigree.
\textit{Berechnen Sie die additiv-genetische Verwandtschaftsmatrix $A$ für das oben angegebene Pedigree}
\points{r lPointsQ1$TaskA}
\end{enumerate}
\solstart
The numerator relationship matrix is computed using pedigreemm::getA(). In a first step, we have to extend the pedigree to contain the founder animals in the ID-column.
vec_founder_sire <- setdiff(tbl_ped_nrm$Sire, tbl_ped_nrm$ID) n_nr_founder_sire <- length(vec_founder_sire) vec_founder_dam <- setdiff(tbl_ped_nrm$Dam, tbl_ped_nrm$ID) n_nr_founder_dam <- length(vec_founder_dam) # check that founder_sire and founder_dam are not the same if (length(intersect(vec_founder_sire, vec_founder_dam)) != 0) stop(" * ERROR: Founder sires and founder dams are not exclusive") tbl_ped_nrm_ext <- dplyr::bind_rows( tibble::tibble(ID = vec_founder_sire[order(vec_founder_sire)], Sex = rep('M', n_nr_founder_sire), Sire = rep(NA, n_nr_founder_sire), Dam = rep(NA, n_nr_founder_sire)), tibble::tibble(ID = vec_founder_dam[order(vec_founder_dam)], Sex = rep('F', n_nr_founder_dam), Sire = rep(NA, n_nr_founder_dam), Dam = rep(NA, n_nr_founder_dam)), tbl_ped_nrm) ped_nrm <- pedigreemm::pedigree(sire = tbl_ped_nrm_ext$Sire, dam = tbl_ped_nrm_ext$Dam, label = as.character(tbl_ped_nrm_ext$ID)) (matA_nrm <- as.matrix(pedigreemm::getA(ped = ped_nrm)))
\solend
\clearpage \pagebreak
\begin{enumerate} \item[b)] Compute the inbreeding coefficients of all animals in the given pedigree. Complete the following table and indicate which of the animals are inbred.
\textit{Berechnen Sie die Inzuchtkoeffizienten aller Tiere im gegebenen Pedigree. Vervollständigen Sie die folgende Tabelle und geben an, welche Tiere ingezüchtet sind.}
\points{r lPointsQ1$TaskB}
\end{enumerate}
\vspace{3ex}
n_nr_animal <- max(tbl_ped_nrm$ID) tbl_inb_empty <- tibble::tibble(ID = 1:n_nr_animal, `Inbreeding Coefficient` = rep('', n_nr_animal), `Animal Inbred (TRUE/FALSE)` = rep('', n_nr_animal)) knitr::kable(tbl_inb_empty, booktabs = TRUE, escape = FALSE, format = 'latex')
\solstart
\vspace{5ex} The numeric solution is
vec_inbr_yn <- ifelse(diag(matA_nrm) > 1, "TRUE", "FALSE") vec_inbr_coef <- diag(matA_nrm) - 1 # tibble for table (tbl_inb_sol <- tibble::tibble(ID = tbl_ped_nrm_ext$ID, `Inbreeding Coefficient` = round(vec_inbr_coef, digits = 3), `Animal Inbred (TRUE/FALSE)` = vec_inbr_yn))
\solend
\clearpage \pagebreak
\begin{enumerate}
\item[c)] Assume that dam r dam_id and sire r sire_id are mated. What is the inbreeding coefficient of their offspring?
\textit{Wir nehmen an, dass die Mutter r dam_id mit dem Vater r sire_id angepaart wird. Wie gross ist der Inzuchtkoeffizient des Nachkommens aus dieser Paarung?}
\points{r lPointsQ1$TaskC}
\end{enumerate}
\solstart
The inbreeding coefficient of the offspring of dam r dam_id and sire r sire_id is half of the relationship between the parents.
0.5 * matA_nrm[sire_id, dam_id]
\solend
\clearpage \pagebreak
tbl_data_pbv <- tibble::tibble(ID = c(4,5,6,7,8), Sire = c(1,3,4,4,6), Dam = c(2,2,5,5,7), Herd = c('Planta','Moos','Moos','Moos','Moos'), `Phenotypic Observation` = c(7.53,8.48,0.26,6.6,2.44)) n_var_u <- 7.2 n_var_e <- 28.8 n_h2 <- 0.2
cat(cnt$out(ps_suffix = "Prediction of Breeding Values"), "\n")
The following dataset is used to predict breeding values.
\textit{Der folgende Datensatz soll für die Schätzung von Zuchtwerten verwendet werden.}
knitr::kable(tbl_data_pbv, booktabs = TRUE, escape = FALSE, format = 'latex')
Assumptions
- the error variance is: $\sigma_e^2 =
r n_var_e$ - the heritability is: $h^2 =
r n_h2$ - the genetic additive variance is: $\sigma_u^2 =
r n_var_u$ - the population mean $\mu$ can be taken as the mean of the given phenotypic observations
\textit{Annahmen}
- \textit{die Restvarianz beträgt:} $\sigma_e^2 =
r n_var_e$ - \textit{die Heritabilität beträgt:} $h^2 =
r n_h2$ - \textit{die genetisch-additive Varianz beträgt:} $\sigma_u^2 =
r n_var_u$ - \textit{das Populationsmittel $\mu$ kann als Mittelwert der phänotypischen Beobachtungen angenommen werden.}
\vspace{3ex} \begin{enumerate} \item[a)] Use the regression method to predict breeding values based on own-performance records for the animals in the table given above.
\textit{Verwenden Sie die Regressionsmethode zur Schätzung der Zuchtwerte basierend auf der Eigenleistung der Tiere in der obigen Tabelle.}
\points{r lPointsQ2$TaskA}
\end{enumerate}
\clearpage \pagebreak
\solstart
According to the regression method, the predicted breeding value ($\hat{u}_i$) for animal $i$ is
$$\hat{u}_i = h^2(y_i - \mu)$$ where $y_i$ is the phenotypic observation of animal $i$, $h^2$ is the heritability and $\mu$ is the population mean.
\vspace{3ex}
vec_y <- tbl_data_pbv$`Phenotypic Observation` n_mu <- mean(vec_y) vec_uhat <- n_h2 * (vec_y - n_mu)
tbl_sol <- tibble::tibble(ID = tbl_data_pbv$ID, PBV = vec_uhat) knitr::kable(tbl_sol, booktabs = TRUE, escape = FALSE, format = 'latex')
\solend
\clearpage \pagebreak
\begin{enumerate} \item[b)] Use a BLUP animal model to predict the breeding values for all animals in the pedigree based on the data given in the table above. Specify the model to predict breeding values, name all model components, compute the expected values and the variance-covariance matrices for all random model components. Insert the information from the above table into the model components where possible. Set up the mixed model equations and compute the solutions for the estimates of fixed effects and for the predicted breeding values.
\textit{Verwenden Sie ein BLUP Tiermodell zur Schätzung der Zuchtwerte aller Tiere im Pedigree basierend auf den Daten in der obigen Tabelle. Spezifizieren Sie das Modell für die Schätzung der Zuchtwerte, benennen Sie alle Modellkomponenten, berechnen Sie die Erwartungswerte und die Varianz-Kovarianz Matrizen aller zufälligen Effekte im Modell. Setzen Sie die verfügbaren Information aus dem Datensatz in die Modellkomponenten ein. Stellen Sie die Mischmodellgleichungen auf und berechnen Sie die Schätzungen der fixen Effekte und der Zuchtwerte.}
\points{r lPointsQ2$TaskB}
\end{enumerate}
\solstart
Model
The BLUP animal model corresponds to the following mixed-effects model
$$y = Xb + Zu + e$$ where $y$ is the vector of observations, $b$ is the vector of fixed effects for the two herds, $u$ is the vector of breeding values for all animals in the pedigree and $e$ is the vector of random residuals. The design matrices $X$ and $Z$ link the fixed effects and the breeding values to the observations, respectively.
Expected Values and Variance-Covariance Matrices
Expected values and variance-covariance matrices of the random components $y$, $u$ and $e$
$$ E \left[\begin{array}{c} y \ u\ e\end{array}\right] = \left[\begin{array}{c} Xb \ 0 \ 0\end{array}\right] \text{, } var \left[\begin{array}{c} y \ u \ e\end{array}\right] = \left[\begin{array}{ccc} ZGZ^T+R & ZG & 0 \ GZ^T & G & 0 \ 0 & 0 & R \end{array}\right] $$
Informations in Model Components
The model vectors are
vec_y <- tbl_data_pbv$`Phenotypic Observation` n_nr_obs <- length(vec_y) n_nr_fix <- nlevels(as.factor(tbl_data_pbv$Herd)) vec_b <- rmdhelp::vecGetVecElem(psBaseElement = 'b', pnVecLen = n_nr_fix, psResult = 'latex') n_nr_animal_pedigree <- length(unique(c(tbl_data_pbv$ID, tbl_data_pbv$Sire, tbl_data_pbv$Dam))) vec_u <- rmdhelp::vecGetVecElem(psBaseElement = 'u', pnVecLen = n_nr_animal_pedigree, psResult = 'latex') vec_e <- rmdhelp::vecGetVecElem(psBaseElement = 'e', pnVecLen = n_nr_obs, psResult = 'latex') cat('$$\n') cat(paste(rmdhelp::bcolumn_vector(pvec = vec_y, ps_name = 'y'), collapse = '\n')) cat('\\text{, }') cat(paste(rmdhelp::bcolumn_vector(pvec = vec_b, ps_name = 'b'), collapse = '\n')) cat('\\text{, }') cat(paste(rmdhelp::bcolumn_vector(pvec = vec_u, ps_name = 'u'), collapse = '\n')) cat('\\text{, }') cat(paste(rmdhelp::bcolumn_vector(pvec = vec_e, ps_name = 'e'), collapse = '\n'), '\n') cat('$$\n')
The model matrices
# matrix X model_mat <- model.matrix(lm(`Phenotypic Observation` ~ Herd, data = tbl_data_pbv)) mat_X <- model_mat[,] dimnames(mat_X) <- NULL mat_X[, 1] <- mat_X[, 1] - mat_X[, 2] # matrix Z mat_Z <- cbind(matrix(0, nrow = n_nr_obs, ncol = (n_nr_animal_pedigree-n_nr_obs)), diag(1, nrow = n_nr_obs, ncol = n_nr_obs)) # output cat('$$\n') cat(paste(rmdhelp::bmatrix(pmat = mat_X, ps_name = 'X'), collapse = '\n'), '\n') cat('\\text{, }') cat(paste(rmdhelp::bmatrix(pmat = mat_Z, ps_name = 'Z'), collapse = '\n'), '\n') cat('$$\n')
Mixed Model Equations
\begin{equation} \left[ \begin{array}{lr} X^TX & X^TZ \ Z^TX & Z^TZ + A^{-1}* \lambda \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] \notag \end{equation}
Solution
# nrm A tbl_ped_pbv_ext <- rmdexam::extend_pedigree(ptbl_ped = tbl_data_pbv[,c('ID', 'Sire', 'Dam')]) ped_pbv <- pedigreemm::pedigree(sire = tbl_ped_pbv_ext$Sire, dam = tbl_ped_pbv_ext$Dam, label = as.character(tbl_ped_pbv_ext$ID)) mat_Ainv_pbv <- as.matrix(pedigreemm::getAInv(ped = ped_pbv)) lambda <- n_var_e / n_var_u # coefficient matrix mat_xtx <- crossprod(mat_X) mat_xtz <- crossprod(mat_X, mat_Z) mat_ztx <- crossprod(mat_Z, mat_X) mat_ztz_ainv_lambda <- crossprod(mat_Z) + lambda * mat_Ainv_pbv mat_coef <- rbind(cbind(mat_xtx,mat_xtz), cbind(mat_ztx,mat_ztz_ainv_lambda)) # rhs mat_xty <- crossprod(mat_X, vec_y) mat_zty <- crossprod(mat_Z, vec_y) mat_rhs <- rbind(mat_xty, mat_zty) # solutions (mat_sol <- solve(mat_coef, mat_rhs))
\solend
\clearpage \pagebreak
tbl_data_gnm <- tibble::tibble(Animal = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15), `SNPLGH` = c(-1,1,1,0,-1,0,-1,1,1,0,1,0,1,-1,1), `SNPFS2` = c(0,1,1,0,0,-1,-1,0,0,-1,-1,0,0,1,0), Observation = c(12.7,46,32.7,19.3,14.8,6.7,2.4,33,29.4,4.9,14.4,19.5,25.6,19.1,25.6)) n_nr_snp <- 2
cat(cnt$out(ps_suffix = "Genomics"), "\n")
Given is the following data set of SNP-Genotyping results.
\textit{Gegeben sind die Genotypisierungsresultate in der nachfolgenden Tabelle.}
knitr::kable(tbl_data_gnm, booktabs = TRUE, escape = FALSE, format = 'latex')
\vspace{3ex} \begin{enumerate} \item[a)] Use a marker effect model to estimate the fixed effects of both markers on the observation. Please specify the fixed-effect model that you use, name all the model components and insert the information from the data into the components where possible.
\textit{Verwenden Sie eine Marker-Effekt Modell zur Schätzung der fixen Effekte der beiden Marker auf die Beobachtung. Bitte spezifizieren Sie das fixe Modell, benennen Sie alle Modellkomponenten und fügen Sie die Daten aus der obigen Tabelle in die Modellkomponenten ein, wo das möglich ist.}
\points{r lPointsQ3$TaskA}
\end{enumerate}
\solstart
Model
The fixed effect model to estimate the marker effects is given by
$$y = Xb + e$$
where $y$ is the vector of observations, $b$ is the vector of fixed effects containing intercept and effect of both markers, $e$ is the vector of random error terms.
Information in Model
The model vectors are
vec_y <- tbl_data_gnm$Observation n_nr_obs <- length(vec_y) vec_b <- rmdhelp::vecGetVecElem(psBaseElement = 'b', pnVecLen = n_nr_snp, psResult = 'latex') vec_e <- rmdhelp::vecGetVecElem(psBaseElement = 'e', pnVecLen = n_nr_obs, psResult = 'latex') cat("$$\n") cat(paste(rmdhelp::bcolumn_vector(pvec = vec_y, ps_name = 'y'), collapse = '\n'), '\n') cat('\\text{, }') cat(paste(rmdhelp::bcolumn_vector(pvec = vec_b, ps_name = 'b'), collapse = '\n'), '\n') cat('\\text{, }') cat(paste(rmdhelp::bcolumn_vector(pvec = vec_e, ps_name = 'e'), collapse = '\n'), '\n') cat("$$\n")
Marker Effects
The marker effects are estimated using lm()
vec_mrk_names <- grep('SNP', colnames(tbl_data_gnm), value = TRUE) lm_mrk_eff <- lm(as.formula(paste0('Observation ~ ', paste0(vec_mrk_names, collapse = ' + '), collapse = '')), data = tbl_data_gnm) summary(lm_mrk_eff)
\solend
\clearpage \pagebreak
\begin{enumerate} \item[b)] Predict the direct genomic breeding values for all animals of the dataset using the marker effects estimated in Task a).
\textit{ Schätzen Sie die direkt-genomischen Zuchtwerte für alle Tiere im Datensatz unter Verwendung der aus Aufgabe a) geschätzten Markereffekte.}
\points{r lPointsQ3$TaskB}
\end{enumerate}
\solstart
Direct genomic breeding values $\hat{g}_i$ for animal $i$ corresponds to the sums of the marker effects, hence
$$\hat{g}_i = x_i^T \cdot \hat{b}$$ where $\hat{b}$ is the vector of estimated marker effects and $x_i^T$ is the row $i$ of the design matrix $X$ corresponding to animal $i$.
vec_mrk_eff <- coefficients(lm_mrk_eff) n_intercept_pos <- which(names(vec_mrk_eff) == '(Intercept)') if (length(n_intercept_pos) > 0) vec_mrk_eff <- vec_mrk_eff[-n_intercept_pos] vec_snp_names <- names(vec_mrk_eff) mat_geno_gnm <- as.matrix(tbl_data_gnm[, vec_snp_names]) (vec_hat_g <- crossprod(t(mat_geno_gnm), vec_mrk_eff))
\solend
\clearpage \pagebreak
n_nr_base_animals <- 280 n_nr_gen <- 16 n_hom_diff <- 50 n_maf <- 0.2 n_level_inb <- 0.01 n_inb_dep <- 0.2 n_dom_dev <- 30
cat(cnt$out(ps_suffix = "Variance and Inbreeding"), "\n")
In dairy cattle semen and embryos of the breeds Brown Swiss and Holstein are often imported from North America. For the Brown Swiss breed, the North American population is based on r n_nr_base_animals female animals. The following assumptions can be made.
- Although, cattle does not follow the same mode of inheritance as the organism shown in the lecture, the computations as shown in the lecture can be used as an approximation.
- In contrast to the size $N$ of the subpopulations that was assumed to be the number of individuals, here $N$ is the number of female animals in a given subpopulation.
\textit{Samen und Embryos der Rassen Brown Swiss und Holstein werden oft aus Nordamerika importiert. Für die Rasse Brown Swiss basiert die Nordamerikanische Population auf r n_nr_base_animals weiblichen Tieren. Die folgenden Annahmen können getroffen werden.}
- \textit{Obwohl das Rind nicht den gleichen Vererbungsmodus zeigt, wie die Organismen, welche in den Vorlesungsunterlagen verwendet wurden, können die dort eingeführten Berechnungen der Inzucht als Annäherungen verwendet werden.}
- \textit{Die Populationsgrösse $N$ hier entspricht der Anzahl weiblichen Tiere.}
\vspace{3ex}
\begin{enumerate}
\item[a)] What is the expected level of inbreeding ($F$) of imported semen in the Brown Swiss breed after r n_nr_gen generations?
\textit{Welches ist der erwartet Wert der Inzucht ($F$) von importiertem Samen in der Rasse Brown Swiss nach r n_nr_gen Generationen?}
\points{r lPointsQ4$TaskA}
\end{enumerate}
\solstart
The level of inbreeding $F_t$ after $t$ generations is computed as
$$F_t = 1 - (1 - \Delta F)^t$$ with $\Delta F = 1/(2N)$ and $N$ being the number of female animals in the base population.
(inb_level_ft <- 1 - (1 - 1/(2*n_nr_base_animals))^n_nr_gen)
\solend
\clearpage \pagebreak
\begin{enumerate}
\item[b)] Compute the between-line, the within-line and the total genetic variance for a single additive Locus where the difference between the homozygous genotypes is r n_hom_diff, the allele frequency $p = r n_maf$ and the level of inbreeding is r n_level_inb.
\textit{Berechnen Sie die Innerhalb-Linie, Zwischen-Linie und die totale genetische Varianz für einen additiven Lokus, wobei die Differenz der homozygoten Genotypen r n_hom_diff entspricht, die Allelefrequenz $p = r n_maf$ ist und der Inzuchtwert r n_level_inb ist.}
\points{r lPointsQ4$TaskB}
\end{enumerate}
\solstart
The variance according to the sources are computed as
tbl_gen_anova <- tibble::tibble(Source = c("Between lines", "Within lines","Total"), Variance = c("$2FV_G$", "$(1-F)V_G$", "$(1+F)V_G$")) knitr::kable(tbl_gen_anova, booktabs = TRUE, escape = FALSE, format = 'latex')
n_vg <- 2 * n_maf * (1-n_maf) * (n_hom_diff/2)^2
with
$$V_G = 2pqa^2 = 2 * r n_maf * r 1-n_maf * r (n_hom_diff/2)^2 = r n_vg$$
Inserting the values, we get
tbl_gen_anova_result <- dplyr::bind_cols(tbl_gen_anova, tibble::tibble(Results = c(2*n_level_inb*n_vg, (1-n_level_inb)*n_vg, (1+n_level_inb)*n_vg))) knitr::kable(tbl_gen_anova_result, booktabs = TRUE, escape = FALSE, format = 'latex')
\solend
\clearpage \pagebreak
\begin{enumerate}
\item[c)] After how many generations is the expected inbreeding depression bigger than r n_inb_dep in a population with $N = r n_nr_base_animals$ animals. The following assumptions can be made
\textit{Nach wie vielen Generationen ist die erwartete Inzuchtdepression grösser als r n_inb_dep in einer Population von $N = r n_nr_base_animals$ Tieren. Die folgenden Annahmen können getroffen werden.}
\points{r lPointsQ4$TaskC}
\end{enumerate}
- single bi-allelic locus
- minor allele frequency $p =
r n_maf$ -
dominance deviation $d =
r n_dom_dev$ -
\textit{einzelner Locus mit zwei Allelen}
- \textit{Frequenz des seltenen Alleles $p =
r n_maf$} - \textit{Dominanzabweichung $d =
r n_dom_dev$}
\solstart
n_inb_f <- n_inb_dep / (2*n_dom_dev*n_maf*(1-n_maf)) n_nr_dgen <- log(1-n_inb_f) / log(1 - 1/(2*n_nr_base_animals))
Inbreeding depression is computed as
$$M_0 - M_F = 2dp(1-p)F = r n_inb_dep$$
hence
$$F = \frac{M_0 - M_F}{2dp(1-p)} = \frac{r n_inb_dep}{2 * r n_dom_dev * r n_maf * (1- r n_maf)} = r n_inb_f$$
From
$$F_t = 1 - (1 - \Delta F)^t$$
we get
$$(1 - \Delta F)^t = 1-F_t$$
and
$$t = \frac{log(1-F_t)}{log(1 - \Delta F)} = r n_nr_dgen$$
Hence after r ceiling(n_nr_dgen) generations the inbreeding depression is bigger than r n_inb_dep
\solend
\clearpage \pagebreak
tbl_data_anova <- tibble::tibble(Animal = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20), Farm = c(3,1,2,3,1,2,3,1,2,1,2,1,2,2,1,1,1,2,3,3), LiveWeight = c(613,621,630,614,629,611,612,614,606,621,621,623,608,603,589,599,610,595,612,616))
cat(cnt$out(ps_suffix = "Variance Components"), "\n")
We are given the following dataset for the trait live weight (LiveWeight) for cattle.
\textit{Der folgende Datensatz umfasst das Merkmal Lebendgewicht (LiveWeight) von Rindern.}
knitr::kable(tbl_data_anova, booktabs = TRUE, escape = FALSE, format = 'latex')
\begin{enumerate} \item[a)] Compute the estimate of the error variance $\sigma_e^2$ from the residuals of the fixed linear model specified below.
\textit{Schätzen Sie die Fehlervarianz $\sigma_e^2$ basierend auf den Residuen des folgenden fixen Modells.}
\points{r lPointsQ5$TaskA}
\end{enumerate}
The fixed linear model that is used is
$$y = Xb + e$$ where $y$ is the vector of all live weight values, $b$ is the vector of the effects caused by the different farms. The fixed linear model is specified in R by
tbl_data_anova$Farm <- as.factor(tbl_data_anova$Farm) lm_lweight <- lm(LiveWeight ~ 0 + Farm, data = tbl_data_anova)
The resulting effects of the farms are
(vec_coef_lweight <- coefficients(lm_lweight))
\solstart
The esimate of the error variance is computed based on the resiudals. The residuals can be obtained by the function residuals() in R.
vec_res <- residuals(lm_lweight) ssq_res <- crossprod(vec_res) (n_est_res_var <- ssq_res / (nrow(tbl_data_anova)-length(vec_coef_lweight)))
The error standard deviation is
(n_est_res_sd <- sqrt(n_est_res_var))
\solend
\clearpage \pagebreak
\begin{enumerate}
\item[b)] Verify your result from task a) with the output of the summary()-function applied to the result of the lm()-function
\textit{Verifizieren Sie das Resultat aus Aufgabe a) anhand des Outputs der summary-Funktion angewendet auf das Resultat der lm()-Funktion}
\points{r lPointsQ5$TaskB}
\end{enumerate}
\solstart
From the task, we have the result object of the lm()-function which is called lm_lweight. Applying the summary()-method leads to
summary(lm_lweight)
The number next to Residual standard error: corresponds to the estimated value of the error standard deviation.
\solend
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.