In charlotte-ngs/ZLHS2016:
knitr::opts_chunk$set(echo = FALSE, results = 'asis')
devtools::load_all()
r6obj_docstat <- rmddochelper::R6ClassDocuStatus$new()
r6obj_docstat$set_current_status(psVersion = "0.0.901",
psStatus = "Initialisation",
psProject = "ZL_HS_2016")
r6obj_docstat$include_doc_stat(psTitle = "## Document Status")
r6ob_abbrtable <- rmddochelper::R6ClassTableAbbrev$new()
### # include table of abbreviations only, if there are any
if (!r6ob_abbrtable$is_empty_abbr())
r6ob_abbrtable$include_abbr_table(psAbbrTitle = "## Abbreviations")
Aufgabe 1: LDL-Zerlegung der Verwandtschaftsmatrix
Gegeben ist das folgende Pedigree
suppressPackageStartupMessages(require(pedigreemm))
nNrAni <- 5
ped <- pedigree(sire = c(NA,NA,1,1,4),dam = c(NA,NA,2,NA,2), label = as.character(1:nNrAni))
print(ped)
a) Stellen Sie die Verwandtschaftsmatrix $A$ auf
b) Ermitteln Sid die LDL-Zerlegung von $A$
c) Berechnen Sie die Inverse der Verwandtschaftsmatrix
Hinweise:
-
Für Pedigrees ohne Inzucht, können die Diagonalelemente $d_{ii}$ nur drei verschiedene Werte annehmen. Diese sind entweder $1$, falls das Tier ein Foundertier ist ohne bekannte Eltern, $3/4$ falls ein Elternteil bekannt ist, oder $1/2$ falls beide Eltern bekannt sind.
-
Ohne Berücksichtigung der Inzucht kann die Inverse $A^{-1}$ einfach aufgrund folgender Regeln aufgestellt werden.
- Initialisiere alle Elemente der Matrix $A^{-1}$ mit $0$
- Für Tier $i$ mit Eltern $s$ und $d$,
- addiere $\delta_i$ zum Element $(i,i)$,
- addiere $-\delta_i/2$ zu den Elementen $(s,i)$, $(i,s)$, $(d,i)$ und $(i,d)$ und
- addiere $\delta_i/4$ zu den Elementen $(s,s)$, $(s,d)$, $(d,s)$ und $(d,d)$
- Für Tier $i$ mit bekanntem Elternteil $d$,
- addiere $\delta_i$ zum Element $(i,i)$,
- addiere $-\delta_i/2$ zu den Elementen $(d,i)$ und $(i,d)$ und
- addiere $\delta_i/4$ zu den Elementen $(d,d)$
- Für Tier $i$ mit unbekannten Eltern
- addiere $\delta_i$ zum Element $(i,i)$
wobei $\delta_i$ das $i$-te Element auf der Diagonalen von $D^{-1}$ ist.
Lösung:
a) Die Verwantschaftsmatrix
matA <- as.matrix(getA(ped=ped))
matR <- t(chol(matA))
cat("$$A = \\left[")
cat(paste(sGetTexMatrix(pmatAMatrix = matA, pnDigits = 4), collapse = "\n"))
cat("\\right]\n")
cat("$$\n")
b) Die LDL-Zerlegung von $A$: Da es im Pedigree keine ingezüchteten Tiere gibt, ist die Matrix $D$ einfach zu bestimmen als.
matD <- diag(Dmat(ped = ped))
cat("$$D = \\left[")
cat(paste(sGetTexMatrix(pmatAMatrix = matD, pnDigits = 4), collapse = "\n"))
cat("\\right]\n")
cat("$$\n")
Die Diagonalelemente sind $1$ für Foundertiere, $3/4$ für Tiere mit einem bekannten Elternteil und $1/2$ für Tiere mit bekannten Eltern.
Die Matrix $L$ ist eine linke untere Dreiecksmatrix, wobei Zeile $i$ dem Mittelwert der Zeilen $s$ und $d$ entsprechen, angenommen, dass die Tiere $s$ und $d$ die Eltern von Tier $i$ sind.
matSInv <- solve(sqrt(matD))
matL <- matR %*% matSInv
cat("$$L = \\left[")
cat(paste(sGetTexMatrix(pmatAMatrix = matL, pnDigits = 4), collapse = "\n"))
cat("\\right]\n")
cat("$$\n")
Daraus finden wir
matDInv <- solve(matD)
cat("$$D^{-1} = \\left[")
cat(paste(sGetTexMatrix(pmatAMatrix = matDInv, pnDigits = 4), collapse = "\n"))
cat("\\right]\n")
cat("$$\n")
und
matLInv <- solve(matL)
cat("$$L^{-1} = \\left[")
cat(paste(sGetTexMatrix(pmatAMatrix = matLInv, pnDigits = 4), collapse = "\n"))
cat("\\right]\n")
cat("$$\n")
Multiplizieren wir die Inversen, dann folgt
matAInv <- as.matrix(getAInv(ped = ped))
cat("$$A^{-1} = (L^{-1})^T * D^{-1} * L^{-1}= \\left[")
cat(paste(sGetTexMatrix(pmatAMatrix = matAInv, pnDigits = 4), collapse = "\n"))
cat("\\right]\n")
cat("$$\n")
\pagebreak
Aufgabe 2: Zuchtwert aufgrund Nachkommenleistung
### # Annahmen
nNrOffspring <- 10
nMeanWwg <- 1.25
nSdWwg <- 0.3
nPopMean <- 1
nH2 <- 0.35
k <- (4-nH2)/nH2
nPhenSd <- 1.1
### # Daten generieren
set.seed(123)
dfWeanWeightGain = data.frame(Tier=1:nNrOffspring, Zuwachs = round(rnorm(nNrOffspring, mean=nMeanWwg, sd=nSdWwg), digits = 2))
nMeanZuw <- mean(dfWeanWeightGain$Zuwachs)
### # ComputeHatA
nHatA <- 2*nNrOffspring/(nNrOffspring + k) * (nMeanZuw-nPopMean)
### # CompSigmaHatA
nSigmaHatA <- (2*nNrOffspring)/(nNrOffspring + k) * sqrt(0.25*nH2 + (1-0.25*nH2)/nNrOffspring) * nPhenSd
### # Obere Grenze
nPvalHatA <- 0.025
nOgHatA <- nHatA + qnorm(nPvalHatA, lower.tail = FALSE) * nSigmaHatA
Stier Elvis hat $r nNrOffspring$ Nachkommen. Von diesen Nachkommen wurden der tägliche Zuwachs (in kg) bis zum Absetzen aufgezeichnet. Das Merkmal hat ein Populationsmittel von $r nPopMean$ kg pro Tag. Die Heritabilität $h^2$ beträgt $r nH2$ und die phänotypische Standardabweichung $\sigma_y$ ist r nPhenSd kg pro Tag.
knitr::kable(dfWeanWeightGain)
a) Schätzen Sie den Zuchtwert von Elvis für das Merkmal tägliche Zunahme bis zum Absetzen aufgrund der Nachkommenleistung
b) Wie gross ist das Bestimmtheitsmass des unter a) geschätzten Zuchtwertes
c) Berechnen Sie aufgrund der Standardabweichung des geschätzten Zuchtwertes die Wahrscheinlichkeit, dass der Zuchtwert von Elvis grösser oder gleich $+r round(nOgHatA, digits = 2)$ kg pro Tag ist.
Lösung:
a) Der geschätzte Zuchtwert für Elvis beträgt
$$\hat{a} = \frac{2n}{n+k}(\tilde{y} - \mu)$$
wobei $n$: Anzahl Beobachtungen, $\tilde{y}$ der Nachkommendurchschnitt, $\mu$ das Populationsmittel und $k=\frac{4-h^2}{h^2}$. Setzen wir diese Werte ein, dann folgt
$$\hat{a} = \frac{2*r nNrOffspring}{r nNrOffspring+r k}(r nMeanZuw - r nPopMean)
= r round(nHatA, digits = 2)$$
b) Das Bestimmtheitsmass für den Zuchtwert von Elvis ist
$$B = r_{a,\tilde{y}}^2 = \frac{n}{n+k} = \frac{r nNrOffspring}{r nNrOffspring+r k}
= r round(nNrOffspring/(nNrOffspring+k), digits=3)$$
\pagebreak
c) Für die Standardabweichung des geschätzten Zuchtwertes, berechnen wir zuerst die Varianz $var(\hat{a})$ und ziehen dann die Wurzel.
\begin{eqnarray}
\sqrt{var(\hat{a})} &=& \sqrt{var(b * (\tilde{y} - \mu))} \nonumber\
&=& b * \sqrt{var(\tilde{y})} \nonumber\
&=& \frac{2n}{n+k}\sqrt{({1\over 4}h^2 + (1-{1\over 4}h^2)/n) \sigma_y^2} \nonumber\
&=& \frac{2 * r nNrOffspring}{r nNrOffspring + r k}
* \sqrt{{1\over 4} * r nH2 + (1 - {1\over 4} * r nH2) / r nNrOffspring} * r nPhenSd\nonumber\
&=& r round(nSigmaHatA, digits=3) \nonumber
\end{eqnarray}
Da $r round(nOgHatA, digits = 2) \approx r round(nHatA, digits=2) + r round(qnorm(0.025, lower.tail = FALSE), digits = 2) * r round(nSigmaHatA, digits=3)$ ist, folgt, dass die Wahrscheinlichkeit rund $r nPvalHatA$ oder $r 100*nPvalHatA\%$ beträgt.
Aufgabe 3: Bedingungen und Schleifen (Loops) in R
Schleifen (Loops) erlauben es uns gewissen Statements wiederholt ausführen zu lassen. Will man als Beispiel alle natürlichen Zahlen zwischen $r nLowerLimit$ und $r nUpperLimit$ ausgeben, dann kann das mit folgendem Loop gemacht werden.
nLowerLimit <- 1
nUpperLimit <- 10
for (nIdx in (nLowerLimit:nUpperLimit)){
cat(nIdx,"\n")
}
Die sogenannten if-Bedingungen können verwendet werden um den Programmablauf zu steuern. Sollen zum Beispiel in der Schleife, welche oben gezeigt wurde, nur die geraden Zahlen ausgegeben werden kann das mit der folgenden Bedingung machen.
nLowerLimit <- 1
nUpperLimit <- 10
for (nIdx in (nLowerLimit:nUpperLimit)){
if (identical(nIdx %% 2, 0)) {
cat(nIdx," ist gerade\n")
}
}
Ihre Aufgabe
Finden Sie in folgendem Pedigree alle Tiere, die keine Mutter haben. Verwenden Sie dazu einen Loop über alle Tiere und testen Sie mit einer if-Bedingung, ob die Mutter bekannt ist.
print(ped)
Hinweise:
- Am einfachsten beginnen Sie, indem Sie das Pedigree mit folgendem Statement der Variablen ped zuweisen.
library(pedigreemm)
nNrAni <- 5
ped <- pedigree(sire = c(NA,NA,1,1,4),dam = c(NA,NA,2,NA,2), label = as.character(1:nNrAni))
- In unseren Pedigrees werden die unbekannten Eltern mit
NA kodiert. Verwenden Sie die Funktion is.na() für die Überprüfung, ob die Mutter bekannt ist.
- Angenommen, sie weisen das Pedigree der Variablen
ped zu, dann finden sie die Anzahl Tiere im Pedigree mit dem Ausdruck length(ped@label)
- Die Mutter des Tieres $i$ bekommen Sie mit dem Ausdruck
ped@dam[i]
Lösung:
suppressPackageStartupMessages(require(pedigreemm))
nNrAni <- 5
ped <- pedigree(sire = c(NA,NA,1,1,4),dam = c(NA,NA,2,NA,2), label = as.character(1:nNrAni))
for (idx in 1:length(ped@label)){
if (is.na(ped@dam[idx])) {
cat("Tier ", idx, " hat eine unbekannte Mutter\n")
}
}
r6ob_abbrtable$writeToTsvFile()
charlotte-ngs/ZLHS2016 documentation built on May 13, 2019, 3: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 = FALSE, results = 'asis') devtools::load_all()
r6obj_docstat <- rmddochelper::R6ClassDocuStatus$new() r6obj_docstat$set_current_status(psVersion = "0.0.901", psStatus = "Initialisation", psProject = "ZL_HS_2016") r6obj_docstat$include_doc_stat(psTitle = "## Document Status")
r6ob_abbrtable <- rmddochelper::R6ClassTableAbbrev$new() ### # include table of abbreviations only, if there are any if (!r6ob_abbrtable$is_empty_abbr()) r6ob_abbrtable$include_abbr_table(psAbbrTitle = "## Abbreviations")
Aufgabe 1: LDL-Zerlegung der Verwandtschaftsmatrix
Gegeben ist das folgende Pedigree
suppressPackageStartupMessages(require(pedigreemm)) nNrAni <- 5 ped <- pedigree(sire = c(NA,NA,1,1,4),dam = c(NA,NA,2,NA,2), label = as.character(1:nNrAni)) print(ped)
a) Stellen Sie die Verwandtschaftsmatrix $A$ auf b) Ermitteln Sid die LDL-Zerlegung von $A$ c) Berechnen Sie die Inverse der Verwandtschaftsmatrix
Hinweise:
-
Für Pedigrees ohne Inzucht, können die Diagonalelemente $d_{ii}$ nur drei verschiedene Werte annehmen. Diese sind entweder $1$, falls das Tier ein Foundertier ist ohne bekannte Eltern, $3/4$ falls ein Elternteil bekannt ist, oder $1/2$ falls beide Eltern bekannt sind.
-
Ohne Berücksichtigung der Inzucht kann die Inverse $A^{-1}$ einfach aufgrund folgender Regeln aufgestellt werden.
- Initialisiere alle Elemente der Matrix $A^{-1}$ mit $0$
- Für Tier $i$ mit Eltern $s$ und $d$,
- addiere $\delta_i$ zum Element $(i,i)$,
- addiere $-\delta_i/2$ zu den Elementen $(s,i)$, $(i,s)$, $(d,i)$ und $(i,d)$ und
- addiere $\delta_i/4$ zu den Elementen $(s,s)$, $(s,d)$, $(d,s)$ und $(d,d)$
- Für Tier $i$ mit bekanntem Elternteil $d$,
- addiere $\delta_i$ zum Element $(i,i)$,
- addiere $-\delta_i/2$ zu den Elementen $(d,i)$ und $(i,d)$ und
- addiere $\delta_i/4$ zu den Elementen $(d,d)$
- Für Tier $i$ mit unbekannten Eltern
- addiere $\delta_i$ zum Element $(i,i)$
wobei $\delta_i$ das $i$-te Element auf der Diagonalen von $D^{-1}$ ist.
Lösung:
a) Die Verwantschaftsmatrix
matA <- as.matrix(getA(ped=ped)) matR <- t(chol(matA)) cat("$$A = \\left[") cat(paste(sGetTexMatrix(pmatAMatrix = matA, pnDigits = 4), collapse = "\n")) cat("\\right]\n") cat("$$\n")
b) Die LDL-Zerlegung von $A$: Da es im Pedigree keine ingezüchteten Tiere gibt, ist die Matrix $D$ einfach zu bestimmen als.
matD <- diag(Dmat(ped = ped)) cat("$$D = \\left[") cat(paste(sGetTexMatrix(pmatAMatrix = matD, pnDigits = 4), collapse = "\n")) cat("\\right]\n") cat("$$\n")
Die Diagonalelemente sind $1$ für Foundertiere, $3/4$ für Tiere mit einem bekannten Elternteil und $1/2$ für Tiere mit bekannten Eltern.
Die Matrix $L$ ist eine linke untere Dreiecksmatrix, wobei Zeile $i$ dem Mittelwert der Zeilen $s$ und $d$ entsprechen, angenommen, dass die Tiere $s$ und $d$ die Eltern von Tier $i$ sind.
matSInv <- solve(sqrt(matD)) matL <- matR %*% matSInv cat("$$L = \\left[") cat(paste(sGetTexMatrix(pmatAMatrix = matL, pnDigits = 4), collapse = "\n")) cat("\\right]\n") cat("$$\n")
Daraus finden wir
matDInv <- solve(matD) cat("$$D^{-1} = \\left[") cat(paste(sGetTexMatrix(pmatAMatrix = matDInv, pnDigits = 4), collapse = "\n")) cat("\\right]\n") cat("$$\n")
und
matLInv <- solve(matL) cat("$$L^{-1} = \\left[") cat(paste(sGetTexMatrix(pmatAMatrix = matLInv, pnDigits = 4), collapse = "\n")) cat("\\right]\n") cat("$$\n")
Multiplizieren wir die Inversen, dann folgt
matAInv <- as.matrix(getAInv(ped = ped)) cat("$$A^{-1} = (L^{-1})^T * D^{-1} * L^{-1}= \\left[") cat(paste(sGetTexMatrix(pmatAMatrix = matAInv, pnDigits = 4), collapse = "\n")) cat("\\right]\n") cat("$$\n")
\pagebreak
Aufgabe 2: Zuchtwert aufgrund Nachkommenleistung
### # Annahmen nNrOffspring <- 10 nMeanWwg <- 1.25 nSdWwg <- 0.3 nPopMean <- 1 nH2 <- 0.35 k <- (4-nH2)/nH2 nPhenSd <- 1.1 ### # Daten generieren set.seed(123) dfWeanWeightGain = data.frame(Tier=1:nNrOffspring, Zuwachs = round(rnorm(nNrOffspring, mean=nMeanWwg, sd=nSdWwg), digits = 2)) nMeanZuw <- mean(dfWeanWeightGain$Zuwachs) ### # ComputeHatA nHatA <- 2*nNrOffspring/(nNrOffspring + k) * (nMeanZuw-nPopMean) ### # CompSigmaHatA nSigmaHatA <- (2*nNrOffspring)/(nNrOffspring + k) * sqrt(0.25*nH2 + (1-0.25*nH2)/nNrOffspring) * nPhenSd ### # Obere Grenze nPvalHatA <- 0.025 nOgHatA <- nHatA + qnorm(nPvalHatA, lower.tail = FALSE) * nSigmaHatA
Stier Elvis hat $r nNrOffspring$ Nachkommen. Von diesen Nachkommen wurden der tägliche Zuwachs (in kg) bis zum Absetzen aufgezeichnet. Das Merkmal hat ein Populationsmittel von $r nPopMean$ kg pro Tag. Die Heritabilität $h^2$ beträgt $r nH2$ und die phänotypische Standardabweichung $\sigma_y$ ist r nPhenSd kg pro Tag.
knitr::kable(dfWeanWeightGain)
a) Schätzen Sie den Zuchtwert von Elvis für das Merkmal tägliche Zunahme bis zum Absetzen aufgrund der Nachkommenleistung
b) Wie gross ist das Bestimmtheitsmass des unter a) geschätzten Zuchtwertes
c) Berechnen Sie aufgrund der Standardabweichung des geschätzten Zuchtwertes die Wahrscheinlichkeit, dass der Zuchtwert von Elvis grösser oder gleich $+r round(nOgHatA, digits = 2)$ kg pro Tag ist.
Lösung:
a) Der geschätzte Zuchtwert für Elvis beträgt
$$\hat{a} = \frac{2n}{n+k}(\tilde{y} - \mu)$$
wobei $n$: Anzahl Beobachtungen, $\tilde{y}$ der Nachkommendurchschnitt, $\mu$ das Populationsmittel und $k=\frac{4-h^2}{h^2}$. Setzen wir diese Werte ein, dann folgt
$$\hat{a} = \frac{2*r nNrOffspring}{r nNrOffspring+r k}(r nMeanZuw - r nPopMean)
= r round(nHatA, digits = 2)$$
b) Das Bestimmtheitsmass für den Zuchtwert von Elvis ist
$$B = r_{a,\tilde{y}}^2 = \frac{n}{n+k} = \frac{r nNrOffspring}{r nNrOffspring+r k}
= r round(nNrOffspring/(nNrOffspring+k), digits=3)$$
\pagebreak
c) Für die Standardabweichung des geschätzten Zuchtwertes, berechnen wir zuerst die Varianz $var(\hat{a})$ und ziehen dann die Wurzel.
\begin{eqnarray}
\sqrt{var(\hat{a})} &=& \sqrt{var(b * (\tilde{y} - \mu))} \nonumber\
&=& b * \sqrt{var(\tilde{y})} \nonumber\
&=& \frac{2n}{n+k}\sqrt{({1\over 4}h^2 + (1-{1\over 4}h^2)/n) \sigma_y^2} \nonumber\
&=& \frac{2 * r nNrOffspring}{r nNrOffspring + r k}
* \sqrt{{1\over 4} * r nH2 + (1 - {1\over 4} * r nH2) / r nNrOffspring} * r nPhenSd\nonumber\
&=& r round(nSigmaHatA, digits=3) \nonumber
\end{eqnarray}
Da $r round(nOgHatA, digits = 2) \approx r round(nHatA, digits=2) + r round(qnorm(0.025, lower.tail = FALSE), digits = 2) * r round(nSigmaHatA, digits=3)$ ist, folgt, dass die Wahrscheinlichkeit rund $r nPvalHatA$ oder $r 100*nPvalHatA\%$ beträgt.
Aufgabe 3: Bedingungen und Schleifen (Loops) in R
Schleifen (Loops) erlauben es uns gewissen Statements wiederholt ausführen zu lassen. Will man als Beispiel alle natürlichen Zahlen zwischen $r nLowerLimit$ und $r nUpperLimit$ ausgeben, dann kann das mit folgendem Loop gemacht werden.
nLowerLimit <- 1 nUpperLimit <- 10 for (nIdx in (nLowerLimit:nUpperLimit)){ cat(nIdx,"\n") }
Die sogenannten if-Bedingungen können verwendet werden um den Programmablauf zu steuern. Sollen zum Beispiel in der Schleife, welche oben gezeigt wurde, nur die geraden Zahlen ausgegeben werden kann das mit der folgenden Bedingung machen.
nLowerLimit <- 1 nUpperLimit <- 10 for (nIdx in (nLowerLimit:nUpperLimit)){ if (identical(nIdx %% 2, 0)) { cat(nIdx," ist gerade\n") } }
Ihre Aufgabe
Finden Sie in folgendem Pedigree alle Tiere, die keine Mutter haben. Verwenden Sie dazu einen Loop über alle Tiere und testen Sie mit einer if-Bedingung, ob die Mutter bekannt ist.
print(ped)
Hinweise:
- Am einfachsten beginnen Sie, indem Sie das Pedigree mit folgendem Statement der Variablen ped zuweisen.
library(pedigreemm) nNrAni <- 5 ped <- pedigree(sire = c(NA,NA,1,1,4),dam = c(NA,NA,2,NA,2), label = as.character(1:nNrAni))
- In unseren Pedigrees werden die unbekannten Eltern mit
NAkodiert. Verwenden Sie die Funktionis.na()für die Überprüfung, ob die Mutter bekannt ist. - Angenommen, sie weisen das Pedigree der Variablen
pedzu, dann finden sie die Anzahl Tiere im Pedigree mit dem Ausdrucklength(ped@label) - Die Mutter des Tieres $i$ bekommen Sie mit dem Ausdruck
ped@dam[i]
Lösung:
suppressPackageStartupMessages(require(pedigreemm)) nNrAni <- 5 ped <- pedigree(sire = c(NA,NA,1,1,4),dam = c(NA,NA,2,NA,2), label = as.character(1:nNrAni))
for (idx in 1:length(ped@label)){ if (is.na(ped@dam[idx])) { cat("Tier ", idx, " hat eine unbekannte Mutter\n") } }
r6ob_abbrtable$writeToTsvFile()
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.