In charlotte-ngs/lbgfs2020: Livestock Breeding and Genomics (Fall 2020)
Genetic Covariances Between Animals {#geneticcovariancesbetweenanimals}
s_this_rmd_file <- rmdhelp::get_this_rmd_file()
mrmt$set_this_rmd_file(ps_this_rmd_file = s_this_rmd_file)
The prediction of breeding values using BLUP as shown in chapter \@ref(blup) uses linear mixed effects models where the breeding value of each animal is included as a random effect. Linear mixed effect models in general and specifically Henderson's mixed model equations require us to be able to specify the variance-covariance matrix of all random effects. When using the animal model, the breeding value of each animal is included as a random effect in the linear mixed effects model. As a consequence of that we need to determine the covariance between the true breeding values of all animals. Figure \@ref(fig:animalcov) tries to display the structure of the required variance-covariance diagrammatically.
#rmddochelper::use_odg_graphic(ps_path = "odg/animalcov.odg")
knitr::include_graphics(path = "odg/animalcov.png")
The variance-covariance matrix shown at the bottom of Figure \@ref(fig:animalcov) has the variances of the true breeding values on the diagonal and all the covariances between the true breeding values of all animals as offdiagonal elements. From the formulation of the linear mixed effect model in \@ref(eq:linearmixedeffectmodelmatrixvector), we defined the variance-covariance matrix of the random effects $u$ to be $G$ (see equation \@ref(eq:variancerandomeffects)). When predicting breeding values with the animal model, the variance-covariance matrix of all components in the vector $u$ is defined as
\begin{equation}
var(u) = G = A * \sigma_u^2
(#eq:genvarcovarmatrixanimalmodel)
\end{equation}
where the matrix $A$ is called numerator relationship matrix.
Similarity Between Individuals {#similaritybetweenindividuals}
At the genetic level there are two different kinds of similarity
- Identity by descent (IBD)
- Identity by state
#rmddochelper::use_odg_graphic(ps_path = "odg/ibdibs.odg")
knitr::include_graphics(path = "odg/ibdibs.png")
Figure \@ref(fig:ibdibs) illustrates the difference between the two type of identities. The type of graph shown in Figure \@ref(fig:ibdibs) is called a pedigree which is commonly used to display the relationship between animals in a population. The rectangle symbols denote male animals and the round symbols stand for female animals. The horizontal connections between female and male animals denote a mating. All animals connected to a vertical line and follow below are progeny of the connected parents.
The notations inside of the symbols stand for the different genotypes of the animals on a given locus. The red arrows denote the path of two $A_1$-alleles which are copies of the same ancestral allele. These two copies are called identical by descent (IBD). The green arrows show the path of two alleles which are identical by state which do not originate from the same copy of any given ancestral alleles.
Numerator Relationship Matrix {#numeratorrelationshipmatrix}
The probability of identical genes by descent (IBD) occurring in two individuals is termed the co-ancestry or the coefficient of kinship [r mrmt$add('Falconer1996')]. The additive genetic relationship between two individuals is twice their co-ancestry. The matrix that expresses the additive genetic relationship among animals in a population is called the numerator relationship matrix $A$. The matrix $A$ is symmetric and its diagonal elements $(A){ii}$ are equal to $1 + F_i$ where $F_i$ is the coefficient of inbreeding of animal $i$. The coefficient of inbreeding $F_i$ indicates whether an animal $i$ is inbred or not. $F_i$ is defined to be half the additive genetic relationship between the parents of $i$. Hence the diagonal element $(A){ii}$ of matrix $A$ corresponds to twice the probability that two gametes taken at random from an animal $i$ will carry IBD-alleles.
The off-diagonal elements $(A)_{ij}$ equals the numerator of the coefficient of relationship between animals $i$ and $j$. Multiplying the matrix $A$ by the additive genetic variance $\sigma_u^2$ leads to the covariance among breeding values. Thus if $u_i$ is the breeding value of animal $i$ then
\begin{equation}
var(u_i) = (A)_{ii} \sigma_u^2 = (1 + F_i) \sigma_u^2
(#eq:vartruebreedingvalue)
\end{equation}
Algorithm To Compute $A$ {#algorithmtocomputea}
The matrix $A$ can be computed using either the
- path coefficient method or
- recursive method.
The second method is especially suitable for an implementation by a software program. In what follows the recursive method to compute the components of $A$ is described now. Initially, animals in a pedigree are numbered from $1$ to $n$ and ordered such that parents precede their progeny. The following rules are then used to compute the components of $A$.
-
If both parents $s$ and $d$ of animal $i$ are known then
- the diagonal element $(A){ii}$ corresponds to: $(A){ii} = 1 + F_i = 1 + {1\over 2} (A)_{sd}$ and
- the off-diagonal element $(A){ji}$ is computed as: $(A){ji} = {1\over 2} ((A){js} + (A){jd})$
- because $A$ is symmetric $(A){ji} = (A){ij}$
-
If only one parent $s$ is known and assumed unrelated to the mate
- $(A)_{ii} = 1$
- $(A){ij} = (A){ji} = {1\over 2} ((A)_{js}$
-
If both parents are unknown
- $(A)_{ii} = 1$
- $(A){ij} = (A){ji} = 0$
Numeric Example
suppressPackageStartupMessages( library(pedigreemm) )
n_nr_ani_ped <- 6
n_nr_parent <- 2
tbl_ped <- tibble::tibble(Calf = c((n_nr_parent+1):n_nr_ani_ped),
Sire = c(1, 1, 4, 5),
Dam = c(2, NA, 3, 2))
ped <- pedigree(sire = c(rep(NA, n_nr_parent), tbl_ped$Sire), dam = c(rep(NA, n_nr_parent), tbl_ped$Dam), label = as.character(1:n_nr_ani_ped))
matA <- as.matrix(getA(ped = ped))
matAinv <- as.matrix(getAInv(ped = ped))
We are given the following pedigree and we want to compute the matrix $A$ using the recursive method described in \@ref(algorithmtocomputea).
knitr::kable(tbl_ped,
format = 'latex',
booktabs = TRUE,
longtable = TRUE,
caption = "Example Pedigree To Compute Additive Genetic Relationship Matrix")
The first step of the computations of $A$ are the numbering and the ordering of all the animals. This is already done in the pedigree shown in Table \@ref(tab:tabpedexample). The components of $A$ are computed row-by-row starting with $(A)_{11}$.
\begin{align}
(A){11} &= 1 + F_1 = 1 + 0 = 1 \notag \
(A){12} &= 0 = (A){21} \notag \
(A){13} &= {1\over 2} ((A){11} + (A){12}) = 0.5 = (A){31} \notag \
(A){14} &= {1\over 2} (A){11} = 0.5 = (A){14} \notag \
(A){15} &= {1\over 2} (A){14} + (A){13}) = 0.5 = (A){51} \notag \
(A){16} &= {1\over 2} (A){15} + (A)_{12}) = 0.25 \notag
\end{align}
The same computations are also done for all the other components of the matrix $A$. The final result for the matrix looks as follows
cat("$$\n")
# cat("A = \\left[")
# cat(paste(rmddochelper::sConvertMatrixToLaTexArray(pmatAMatrix = matA, pnDigits = 4), sep = "\n"), "\n")
# cat("\\right]\n")
cat(paste(rmdhelp::bmatrix(pmat = matA, ps_name = 'A'), sep = '\n'), '\n')
cat("$$\n")
As a result, we can see from the components of the above shown matrix $A$ that animals $1$ and $2$ are not related to each other. Furthermore from the diagonal elements of $A$, it follows that animals $5$ and $6$ are inbred while animals $1$ to $4$ are not inbred.
\pagebreak
The Inverse Numerator Relationship Matrix {#inversenumeratorrelationshipmatrix}
The general form of Henderson's mixed model equations as shown in \@ref(eq:mixedmodeleq) depend on four matrices
- Design matrix $X$ for the fixed effects
- Design matrix $Z$ for the random effects
- The inverse variance-covariance matrix $R^{-1}$ for the residuals $e$ and
- The inverse variance-covariance matrix $G^{-1}$ for the random breeding values $a$.
When using the animal model as described in section \@ref(#animalmodel), we have seen in \@ref(eq:genvarcovarmatrixanimalmodel) that $G$ corresponds to the product of the numerator relationship matrix $A$ and the genetic additive variance $\sigma_u^2$. But the mixed model equations depend on the inverse $G^{-1}$ of $G$ which means
\begin{equation}
G^{-1} = (A * \sigma_u^2)^{-1} = A^{-1} * {1\over \sigma_u^2}
(#eq:ginvformme)
\end{equation}
From \@ref(eq:ginvformme) we can see that in order to be able to set up the mixed model equations for an animal model, we need the inverse numerator relationship matrix $A^{-1}$. Because in practical routine predictions of breeding values, we have something in the order of 10 million animals that we predict breeding values for. Hence the matrix $A$ has 10 million rows and 10 million columns. A matrix of that size cannot be inverted explicitly with commonly known methods such as described under https://en.wikipedia.org/wiki/Invertible_matrix. This would have been the end of the BLUP animal model, if not Henderson published in [r mrmt$add('Henderson1976')] a set of simple rules to directly construct the matrix $A^{-1}$ without setting up the numerator relationship matrix.
Structure of $A^{-1}$
When looking at a concrete example of an inverse of a numerator relationship matrix as shown below, we can discover some important facts. Let us assume the following pedigree.
### # first six animals from Goetz p. 83
suppressPackageStartupMessages( library(pedigreemm) )
n_nr_ani_ped <- 5
n_nr_parent <- 3
tbl_ped_simple <- dplyr::data_frame(Calf = c(1:n_nr_ani_ped),
Sire = c(NA, NA, NA, 1, 3),
Dam = c(NA, NA, NA, 2, 2))
ped_simple <- pedigree(sire = tbl_ped_simple$Sire, dam = tbl_ped_simple$Dam, label = as.character(1:n_nr_ani_ped))
matu_simple <- as.matrix(getA(ped = ped_simple))
matAinv_simple <- as.matrix(getAInv(ped = ped_simple))
knitr::kable(tbl_ped_simple,
format = 'latex',
booktabs = TRUE,
longtable = TRUE,
caption = "Pedigree Used To Compute Inverse Numerator Relationship Matrix")
The numerator relationship matrix $A$ for the pedigree shown in Table \@ref(tab:simplepedexamplesetup) is shown in \@ref(eq:numeratorrelationshipmatrix).
cat("\\begin{equation}\n")
# cat("A = \\left[")
# cat(paste(rmddochelper::sConvertMatrixToLaTexArray(pmatAMatrix = matu_simple, pnDigits = 4), sep = "\n"), "\n")
# cat("\\right]\n")
cat(paste(rmdhelp::bmatrix(pmat = matu_simple, ps_name = 'A'), sep = '\n'), '\n')
cat("\\label{eq:numeratorrelationshipmatrix}\n")
cat("\\end{equation}\n")
The matrix $A^{-1}$ shown below corresponds to the inverse of the matrix $A$ from \@ref(eq:numeratorrelationshipmatrix).
cat("$$\n")
# cat("A^{-1} = \\left[")
# cat(paste(rmddochelper::sConvertMatrixToLaTexArray(pmatAMatrix = matAinv_simple, pnDigits = 4), sep = "\n"), "\n")
# cat("\\right]")
cat(paste(rmdhelp::bmatrix(pmat = matAinv_simple, ps_name = 'A^{-1}'), sep = '\n'), '\n')
cat("$$")
From the above shown inverse $A^{-1}$, we recognize that $A^{-1}$ has a simpler structure than $A$ itself. The reason for this is that non-zero elements occur only at matrix elements of $A^{-1}$ corresponding to parents and progeny or to mates. Furthermore off-diagonal elements corresponding to mates are positive and elements corresponding to parents and progeny are negative. These observations were used by Henderson in [@Henderson1976] to come up with the rules described below.
Henderson's Rule To Set Up $A^{-1}$
We denote the row or column index corresponding to an animal of interest as $i$ and the row or column index belonging to the animals father as $s$ and the row or column index corresponding to animal $i$'s mother as $d$. The rules differentiate the following three cases
- both parents of animal $i$ are known
- only one parent of animal $i$ is known
- both parents of animal $i$ are unknown
Both Parents Known
- add $2$ to the diagonal-element $(i,i)$
- add $-1$ to off-diagonal elements $(s,i)$, $(i,s)$, $(d,i)$ and $(i,d)$
- add ${1\over 2}$ to elements $(s,s)$, $(d,d)$, $(s,d)$, $(d,s)$
Only One Parent Known
We assume that sire $s$ is known
- add ${4\over 3}$ to diagonal-element $(i,i)$
- add $-{2\over 3}$ to off-diagonal elements $(s,i)$, $(i,s)$
- add ${1\over 3}$ to element $(s,s)$
Both Parents Unknown
- add $1$ to diagonal-element $(i,i)$
The application of Henderson's rules to construct $A^{-1}$ directly will be a problem in one of the coming exercises. When applying the above described rules, it has to be noted well that this simple version of the rules are only valid for a pedigree without inbreeding. We will see in section \@ref(#derivationofhendersonsrules) how to extend the simple rules such that they can be used for a general pedigree with inbreeding.
Derivation of Henderson's Rules {#derivationofhendersonsrules}
Henderson's rules can be related to the so-called LDL-decomposition of the numerator relationship matrix $A$. The result of this decomposition consists of the equivalence between matrix $A$ and the product of three matrices $L$, $D$ and $L^T$. The matrix $L$ is a lower triangular matrix and the matrix $D$ is a diagonal matrix. The reason for why we are doing this decomposition of $A$ is that the matrices $L$ and $D$ can much easier be inverted than the matrix $A$. The LDL-decomposition is a general procedure in linear algebra and it can be applied to any symmetric and positive-definite matrix not only to numerator relationship matrices. But when the LDL-decomposition is applied to a numerator relationship matrix, the matrices $L$ and $D$ do also have a special genetic meaning. This meaning is revealed in the following derivation.
Decomposition of True Breeding Value and its Variance {#decompositionoftruebreedingvalueanditsvariance}
The true breeding value ($u_i$) of animal $i$ can be decomposed into the true breeding values of animal $i$'s parents $s$ and $d$ and the mendelian sampling term $m_i$
\begin{equation}
u_i = {1\over 2} u_s + {1\over 2} u_d + m_i
(#eq:decomposeanimaltruebreedingvalue)
\end{equation}
Applying the decomposition shown in \@ref(eq:decomposeanimaltruebreedingvalue) for all animals in the pedigree and combining the decompositions into a matrix-vector notation, we get
\begin{equation}
u = P \cdot u + m
(#eq:decomposetruebreedingvaluematvec)
\end{equation}
\begin{tabular}{lll}
where & & \
& $u$ & vector of true breeding values for all animals \
& $P$ & matrix linking animals to their parents \
& $m$ & vector of mendelian sampling terms
\end{tabular}
Equation \@ref(eq:decomposetruebreedingvaluematvec) can be interpreted as regression of the true breeding values onto parental breeding values. In such a regression the random term $m$ is like a residual term. The genetic meaning of $m$ corresponds to the deviation of $u_i$ from the full-sib average of the true breeding values within the nuclear family with parents $s$ and $d$. The term $m$ is called mendelian sampling term. The source of $m$ is the fact that during meiosis, parental alleles are randomly assigned to each progeny. Bulmer [r mrmt$add('Bulmer1971')] has shown that $m_i$ are independent from true breeding values of parents $s$ and $d$. Using this result, we take the variance on both sides of \@ref(eq:decomposeanimaltruebreedingvalue) leading to
\begin{align}
var(u_i) &= var({1\over 2} u_s + {1\over 2} u_d + m_i) \notag \
&= {1\over 4} var(u_s) + {1\over 4} var(u_d) + {1\over 2} cov(u_s, u_d) + var(m_i)
(#eq:vartruebreedingvalue)
\end{align}
From \@ref(eq:genvarcovarmatrixanimalmodel) together with the definition of the numerator relationship matrix $A$, we know that
\begin{align}
var(u_i) &= (1 + F_i) \sigma_u^2 \notag \
var(u_s) &= (1 + F_s) \sigma_u^2 \notag \
var(u_d) &= (1 + F_d) \sigma_u^2 \notag \
cov(u_s, u_d) &= (A)_{sd} \sigma_u^2 = 2F_i \sigma_u^2
(#eq:usedefnumeratorrelationshipmatrix)
\end{align}
Variance of Mendelian Sampling Terms {#variancemendeliansamplingterm}
Inserting the relations from \@ref(eq:usedefnumeratorrelationshipmatrix) into \@ref(eq:vartruebreedingvalue) and solving for $var(m_i)$ gives the following result
\begin{align}
var(m_i) &= var(u_i) - {1\over 4} var(u_s) - {1\over 4} var(u_d) - {1\over 2} cov(u_s, u_d) \notag \
&= (1 + F_i) \sigma_u^2 - {1\over 4}(1 + F_s) \sigma_u^2 - {1\over 4} (1 + F_d) \sigma_u^2 - {1\over 2} * 2 * F_i * \sigma_u^2 \notag \
&= \left({1\over 2} - {1\over 4}(F_s + F_d)\right) \sigma_u^2
(#eq:mendeliansamplingvariance)
\end{align}
In case where only parent $s$ of animal $i$ is known the terms in \@ref(eq:decomposeanimaltruebreedingvalue) and \@ref(eq:mendeliansamplingvariance) change to
\begin{equation}
u_i = {1\over 2} u_s + {1\over 2} u_d + m_i
(#eq:decomposeanimaltruebreedingvalue)
\end{equation}
\begin{align}
u_i &= {1\over 2} u_s + m_i \notag \
var(m_i) &= \left(1 - {1\over 4}(1+ F_s)\right) \sigma_u^2 \notag \
&= \left({3\over 4} - {1\over 4}F_s\right) \sigma_u^2
(#eq:mendeliansamplingvarianceoneparentknown)
\end{align}
When both parents are unknown, we get
\begin{align}
u_i &= m_i \notag \
var(m_i) &= \sigma_u^2
(#eq:mendeliansamplingvariancebothparentunknown)
\end{align}
Decomposition of $A$ {#decompositionofa}
The true breeding values $u_s$ of sire $s$ and $u_d$ of dam $d$ can be decomposed in a similar way as shown for the true breeding value $u_i$ in \@ref(eq:decomposeanimaltruebreedingvalue).
\begin{align}
u_s &= {1\over 2} u_{ss} + {1\over 2} u_{ds} + m_s \notag \
u_d &= {1\over 2} u_{sd} + {1\over 2} u_{dd} + m_d
(#eq:decomposeparenttruebreedingvalue)
\end{align}
\begin{tabular}{lll}
where & & \
& $ss$ & sire of $s$ \
& $ds$ & dam of $s$ \
& $sd$ & sire of $d$ \
& $dd$ & dam of $d$
\end{tabular}
Using \@ref(eq:decomposeparenttruebreedingvalue) together with \@ref(eq:decomposeanimaltruebreedingvalue) leads to the following expression for $u_i$
\begin{align}
u_i &= {1\over 2} u_s + {1\over 2} u_d + m_i \notag \
&= {1\over 2} ({1\over 2} u_{ss} + {1\over 2} u_{ds} + m_s)+ {1\over 2} ({1\over 2} u_{sd} + {1\over 2} u_{dd} + m_d) + m_i \notag \
&= {1\over 4} u_{ss} + {1\over 4} u_{ds} + {1\over 4} u_{sd} + {1\over 4} u_{dd} + {1\over 2} m_s + {1\over 2} m_d + m_i \notag
(#eq:recursivedecomposeanimaltruebreedingvalue)
\end{align}
This type of decomposition can also be done for the grand-parents of animal $i$ and further back until we get to true breeding values of animals with unknown parents. Animals of ancestor generations with unknown parents are called founder population. The process of decomposing true breeding values back through all generations of a pedigree is called recursive decomposition of animal $i$'s true breeding value. Because we know from \@ref(eq:mendeliansamplingvariancebothparentunknown) that the decomposition of true breeding values $u_k$ for an animal $k$ with unknown parents is $u_k = m_k$, the result of the recursive decomposition of $u_i$ is a sum of mendelian sampling terms linking the ancestors of $i$ back to the founder population.
Ordering all animals in a pedigree according to their age and writing the result of the recursive decomposition of all true breeding values in matrix-vector notation leads to
\begin{equation}
u = L \cdot m
(#eq:resultrecursivedecomposition)
\end{equation}
The matrix $L$ is a lower triangular matrix. The row corresponding to animal $i$ is the average of the rows in $L$ corresponding to parents $s$ and $d$ of $i$. The matrix $L$ contains the path of the gene flow from the base population to all animals in the population. From equation \@ref(eq:resultrecursivedecomposition), we are computing the variance of all true breeding values which leads to
\begin{align}
var(u) &= var(L \cdot m) \notag \
&= L \cdot var(m) \cdot L^T \notag \
&= L \cdot (D \sigma_u^2) \cdot L^T \notag \
&= L \cdot D \cdot L^T \sigma_u^2 \notag
&= A \sigma_u^2
(#eq:vararecursivedecomposition)
\end{align}
From \@ref(eq:vararecursivedecomposition), the LDL-decomposition of $A$ follows directly as $A = LDL^T$. The single components $m_i$ are independent of each other. This also means that $cov(m_i, m_j) = 0$ for $i \ne j$. Hence the matrix $D$ is a diagonal matrix. In section \@ref(#variancemendeliansamplingterm) we have seen that $var(m_i)$ always contain $\sigma_u^2$ as a factor. Hence it is reasonable to express $var(m)$ as a product of a diagonal matrix $D$ times the genetic additive variance $\sigma_u^2$. The diagonal elements of matrix $D$ are computed as shown in section \@ref(#variancemendeliansamplingterm).
The mixed model equations require the inverse numerator relationship matrix $A^{-1}$. Thanks to the LDL-decomposition of $A$, we can compute $A^{-1}$ as
\begin{equation}
A^{-1} = (L \cdot D \cdot L^T)^{-1} = (L^T)^{-1} \cdot D^{-1} \cdot L^{-1}
(#eq:ainvldldecomp)
\end{equation}
The inverse $D^{-1}$ is easy to compute, because $D$ is a diagonal matrix. As a consequence of that $D^{-1}$ is also a diagonal matrix where the elements $(D^{-1}){ii}$ correspond to the inverse of the elements of the original matrix $D$, i.e. $(D^{-1}){ii} = 1/(D)_{ii}$. The matrix $L^{-1}$ is also a lower triangular matrix and can easily computed based on the two relations for the vector $m$. Based on \@ref(eq:decomposetruebreedingvaluematvec), we know
\begin{equation}
m = u - Pu = (I-P)u
(#eq:mendeliansamplingiminuspa)
\end{equation}
and from \@ref(eq:resultrecursivedecomposition), we get
\begin{equation}
m = L^{-1}u
(#eq:mendeliansamplinglinversea)
\end{equation}
Setting both expressions for $m$ in \@ref(eq:mendeliansamplingiminuspa) and \@ref(eq:mendeliansamplinglinversea) equal can be used to obtain $L^{-1}$
\begin{equation}
m = (I-P)u = L^{-1}u
(#eq:linverseequation)
\end{equation}
Therefore,
\begin{equation}
L^{-1} = I - P
(#eq:linverseresult)
\end{equation}
General Version of Henderson's Rule
Based on the decomposition of $A^{-1}$ shown in \@ref(eq:ainvldldecomp), the general version of Henderson's rule where inbreeding of animals can be accounted for can be summarized as
- Start with a matrix $A^{-1}$ where all elements are set to $0$.
- Let $d^i$ be the $i$-th diagonal element of $D^{-1}$ for animal $i$, assuming $i$ has parents $s$ and $d$.
- Then add the following contributions to $A^{-1}$
- $d^i$ to the element $(i,i)$
- $-d^i/2$ to the elements $(s,i)$, $(i,s)$, $(d,i)$, $(i,d)$
- $d^i/4$ to the elements $(s,s)$, $(s,d)$, $(d,s)$, $(d,d)$
The contributions to rows and columns corresponding to parents $s$ and $d$ are only added, if they are known. Because the elements $d^i$ are dependent on the inbreeding coefficients $F_s$ and $F_d$ of parents $s$ and $d$, we have to find an efficient way to compute inbreeding coefficients for all animals in a population.
Computing Inbreeding Coefficients {#computinginbreedingcoefficients}
Inbreeding coefficients can be computed using different methods. From all these methods, we are just showing the one method which was described in [r mrmt$add('Quaas1976')]. This method is based on the properties of a second decomposition of the numerator relationship matrix $A$ which is called the cholesky decomposition. The cholesky decomposition of a matrix $A$ corresponds to
\begin{equation}
A = RR^T
(#eq:choleskymata)
\end{equation}
where the matrix $R$ is a lower triangular matrix. At this point we have to be clear about this. In practical routine evaluations, we will not explicitly compute this decomposition, because we do not want to construct $A$ explicitly. We are just using the properties of \@ref(eq:choleskymata) to find an efficient way to compute the diagonal elements $(A){ii}$ of $A$ without constructing the complete matrix $A$. The diagonal elements $(A){ii}$ are important, because they contain the inbreeding coefficients ($F_i$), as $(A){ii} = 1 + F_i$. Based on \@ref(eq:choleskymata), $(A){ii}$ can be computed from the components of $R$ as
\begin{equation}
(A){ii} = \sum{j=1}^i (R)_{ij}^2
(#eq:DiagElemMatA)
\end{equation}
This can be shown with a small $3\times 3$ matrix $A$
nAnzAni <- 3
matA <- rmdhelp::matGetMatElem(psBaseElement = "(A)", pnNrRow = nAnzAni, pnNrCol = nAnzAni)
matR <- rmdhelp::matLowerTri(psBaseElement = "(R)", pnNrRow = nAnzAni, pnNrCol = nAnzAni)
cat("$$")
cat(paste(rmdhelp::bmatrix(pmat = matA), collapse = "\n"))
cat(" \n")
cat(" = ")
cat(paste(rmdhelp::bmatrix(pmat = matR), collapse = "\n"))
cat(" \n")
cat(" * ")
cat(paste(rmdhelp::bmatrix(pmat = t(matR)), collapse = "\n"))
cat("\n")
cat("$$\n")
For the above shown example, the diagonal elements $(A)_{ii}$ are computed as
\begin{align}
(A){11} &= (R){11}^2 \notag \
(A){22} &= (R){21}^2 + (R){22}^2 \notag \
(A){33} &= (R){31}^2 + (R){32}^2 + (R)_{33}^2 \notag
(#eq:MatADiag)
\end{align}
This shows that diagonal elements $(A)_{ii}$ can be computed using just all the components of row $i$ in $R$ up to the diagonal. Next, we have to show how to compute the components of $R$.
Recursive Computation of $R$
We are setting the two decompositions of $A$ equal which leads to
\begin{equation}
A = R * R^T = L * D * L^T
(#eq:CholEqLdl)
\end{equation}
Let us write the matrix $R$ as a product of two matrices $L$ and $S$ where $L$ corresponds to the matrix from the LDL-decomposition and insert that product into \@ref(eq:CholEqLdl)
\begin{equation}
A = R * R^T = L * S * (L * S)^T = L * S * S^T * L^T= L * D * L^T
(#eq:CholRLSEqLdl)
\end{equation}
From \@ref(eq:CholRLSEqLdl), we conclude that $D = S \cdot S^T$ where $S$ is also a diagonal matrix with elements $(S){ii} = \sqrt{(D){ii}}$. For our small example we get
\begin{equation}
matL <- rmdhelp::matLowerTri(psBaseElement = "(L)", pnNrRow = nAnzAni, pnNrCol = nAnzAni, pvecDiag = 1)
matS <- rmdhelp::matDiag(psBaseElement = "(S)", pnNrRow = nAnzAni, pnNrCol = nAnzAni)
cat(paste(rmdhelp::bmatrix(pmat = matR), collapse = "\n"))
cat(" = \n")
cat(paste(rmdhelp::bmatrix(pmat = matL), collapse = "\n"))
cat(" * \n")
cat(paste(rmdhelp::bmatrix(pmat = matS), collapse = "\n"))
cat("\n")
(#eq:SmallExRLS)
\end{equation}
From \@ref(eq:SmallExRLS), we can see that the diagonal elements $(R){ii}$ are equal to $(S){ii}$. Therefore
\begin{equation}
(R){ii} = (S){ii} = \sqrt{(D)_{ii}}
\label{eq:ComputeRii}
\end{equation}
Earlier, we have seen that diagonal elements $(D)_{ii}$ of $D$ correspond to
\begin{equation}
(D){ii} = {1\over 2}\ - {1\over 4}\ (F_s + F_d) = 1 - 0.25((A){ss} + (A)_{dd})
\label{eq:ComputeDii}
\end{equation}
and hence
\begin{equation}
(R){ii} = \sqrt{1 - 0.25((A){ss} + (A)_{dd})}
\label{eq:ComputeRiiResult}
\end{equation}
The components $(A){ss}$ and $(A){dd}$ correspond to diagonal elements of $A$ of parents of $s$ and $d$ which were computed earlier.
The off-diagonal elements of $R$ are computed as
\begin{equation}
(R){ij} = (L){ij} * (S)_{jj}
\label{eq:ComputeRij}
\end{equation}
One property of the matrix $L$ is that any element $(L){ij}$ is equal to the average of elements $(L){sj}$ and $(L)_{dj}$, if $s$ and $d$ are parents of animal $i$. Using this we get
\begin{align}
(R){ij} &= (L){ij} * (S){jj} \nonumber\
&= {1\over 2}((L){sj} + (L){dj}) * (S){jj} \nonumber\
&= {1\over 2}((R){sj} + (R){dj})
(#eq:ComputeRijRec)
\end{align}
With that we have shown how to compute all elements of $R$ recursively. This requires an ordering of all animals according to their age.
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.
Genetic Covariances Between Animals {#geneticcovariancesbetweenanimals}
s_this_rmd_file <- rmdhelp::get_this_rmd_file() mrmt$set_this_rmd_file(ps_this_rmd_file = s_this_rmd_file)
The prediction of breeding values using BLUP as shown in chapter \@ref(blup) uses linear mixed effects models where the breeding value of each animal is included as a random effect. Linear mixed effect models in general and specifically Henderson's mixed model equations require us to be able to specify the variance-covariance matrix of all random effects. When using the animal model, the breeding value of each animal is included as a random effect in the linear mixed effects model. As a consequence of that we need to determine the covariance between the true breeding values of all animals. Figure \@ref(fig:animalcov) tries to display the structure of the required variance-covariance diagrammatically.
#rmddochelper::use_odg_graphic(ps_path = "odg/animalcov.odg") knitr::include_graphics(path = "odg/animalcov.png")
The variance-covariance matrix shown at the bottom of Figure \@ref(fig:animalcov) has the variances of the true breeding values on the diagonal and all the covariances between the true breeding values of all animals as offdiagonal elements. From the formulation of the linear mixed effect model in \@ref(eq:linearmixedeffectmodelmatrixvector), we defined the variance-covariance matrix of the random effects $u$ to be $G$ (see equation \@ref(eq:variancerandomeffects)). When predicting breeding values with the animal model, the variance-covariance matrix of all components in the vector $u$ is defined as
\begin{equation} var(u) = G = A * \sigma_u^2 (#eq:genvarcovarmatrixanimalmodel) \end{equation}
where the matrix $A$ is called numerator relationship matrix.
Similarity Between Individuals {#similaritybetweenindividuals}
At the genetic level there are two different kinds of similarity
- Identity by descent (IBD)
- Identity by state
#rmddochelper::use_odg_graphic(ps_path = "odg/ibdibs.odg") knitr::include_graphics(path = "odg/ibdibs.png")
Figure \@ref(fig:ibdibs) illustrates the difference between the two type of identities. The type of graph shown in Figure \@ref(fig:ibdibs) is called a pedigree which is commonly used to display the relationship between animals in a population. The rectangle symbols denote male animals and the round symbols stand for female animals. The horizontal connections between female and male animals denote a mating. All animals connected to a vertical line and follow below are progeny of the connected parents.
The notations inside of the symbols stand for the different genotypes of the animals on a given locus. The red arrows denote the path of two $A_1$-alleles which are copies of the same ancestral allele. These two copies are called identical by descent (IBD). The green arrows show the path of two alleles which are identical by state which do not originate from the same copy of any given ancestral alleles.
Numerator Relationship Matrix {#numeratorrelationshipmatrix}
The probability of identical genes by descent (IBD) occurring in two individuals is termed the co-ancestry or the coefficient of kinship [r mrmt$add('Falconer1996')]. The additive genetic relationship between two individuals is twice their co-ancestry. The matrix that expresses the additive genetic relationship among animals in a population is called the numerator relationship matrix $A$. The matrix $A$ is symmetric and its diagonal elements $(A){ii}$ are equal to $1 + F_i$ where $F_i$ is the coefficient of inbreeding of animal $i$. The coefficient of inbreeding $F_i$ indicates whether an animal $i$ is inbred or not. $F_i$ is defined to be half the additive genetic relationship between the parents of $i$. Hence the diagonal element $(A){ii}$ of matrix $A$ corresponds to twice the probability that two gametes taken at random from an animal $i$ will carry IBD-alleles.
The off-diagonal elements $(A)_{ij}$ equals the numerator of the coefficient of relationship between animals $i$ and $j$. Multiplying the matrix $A$ by the additive genetic variance $\sigma_u^2$ leads to the covariance among breeding values. Thus if $u_i$ is the breeding value of animal $i$ then
\begin{equation} var(u_i) = (A)_{ii} \sigma_u^2 = (1 + F_i) \sigma_u^2 (#eq:vartruebreedingvalue) \end{equation}
Algorithm To Compute $A$ {#algorithmtocomputea}
The matrix $A$ can be computed using either the
- path coefficient method or
- recursive method.
The second method is especially suitable for an implementation by a software program. In what follows the recursive method to compute the components of $A$ is described now. Initially, animals in a pedigree are numbered from $1$ to $n$ and ordered such that parents precede their progeny. The following rules are then used to compute the components of $A$.
-
If both parents $s$ and $d$ of animal $i$ are known then
- the diagonal element $(A){ii}$ corresponds to: $(A){ii} = 1 + F_i = 1 + {1\over 2} (A)_{sd}$ and
- the off-diagonal element $(A){ji}$ is computed as: $(A){ji} = {1\over 2} ((A){js} + (A){jd})$
- because $A$ is symmetric $(A){ji} = (A){ij}$
-
If only one parent $s$ is known and assumed unrelated to the mate
- $(A)_{ii} = 1$
- $(A){ij} = (A){ji} = {1\over 2} ((A)_{js}$
-
If both parents are unknown
- $(A)_{ii} = 1$
- $(A){ij} = (A){ji} = 0$
Numeric Example
suppressPackageStartupMessages( library(pedigreemm) ) n_nr_ani_ped <- 6 n_nr_parent <- 2 tbl_ped <- tibble::tibble(Calf = c((n_nr_parent+1):n_nr_ani_ped), Sire = c(1, 1, 4, 5), Dam = c(2, NA, 3, 2)) ped <- pedigree(sire = c(rep(NA, n_nr_parent), tbl_ped$Sire), dam = c(rep(NA, n_nr_parent), tbl_ped$Dam), label = as.character(1:n_nr_ani_ped)) matA <- as.matrix(getA(ped = ped)) matAinv <- as.matrix(getAInv(ped = ped))
We are given the following pedigree and we want to compute the matrix $A$ using the recursive method described in \@ref(algorithmtocomputea).
knitr::kable(tbl_ped, format = 'latex', booktabs = TRUE, longtable = TRUE, caption = "Example Pedigree To Compute Additive Genetic Relationship Matrix")
The first step of the computations of $A$ are the numbering and the ordering of all the animals. This is already done in the pedigree shown in Table \@ref(tab:tabpedexample). The components of $A$ are computed row-by-row starting with $(A)_{11}$.
\begin{align} (A){11} &= 1 + F_1 = 1 + 0 = 1 \notag \ (A){12} &= 0 = (A){21} \notag \ (A){13} &= {1\over 2} ((A){11} + (A){12}) = 0.5 = (A){31} \notag \ (A){14} &= {1\over 2} (A){11} = 0.5 = (A){14} \notag \ (A){15} &= {1\over 2} (A){14} + (A){13}) = 0.5 = (A){51} \notag \ (A){16} &= {1\over 2} (A){15} + (A)_{12}) = 0.25 \notag \end{align}
The same computations are also done for all the other components of the matrix $A$. The final result for the matrix looks as follows
cat("$$\n") # cat("A = \\left[") # cat(paste(rmddochelper::sConvertMatrixToLaTexArray(pmatAMatrix = matA, pnDigits = 4), sep = "\n"), "\n") # cat("\\right]\n") cat(paste(rmdhelp::bmatrix(pmat = matA, ps_name = 'A'), sep = '\n'), '\n') cat("$$\n")
As a result, we can see from the components of the above shown matrix $A$ that animals $1$ and $2$ are not related to each other. Furthermore from the diagonal elements of $A$, it follows that animals $5$ and $6$ are inbred while animals $1$ to $4$ are not inbred.
\pagebreak
The Inverse Numerator Relationship Matrix {#inversenumeratorrelationshipmatrix}
The general form of Henderson's mixed model equations as shown in \@ref(eq:mixedmodeleq) depend on four matrices
- Design matrix $X$ for the fixed effects
- Design matrix $Z$ for the random effects
- The inverse variance-covariance matrix $R^{-1}$ for the residuals $e$ and
- The inverse variance-covariance matrix $G^{-1}$ for the random breeding values $a$.
When using the animal model as described in section \@ref(#animalmodel), we have seen in \@ref(eq:genvarcovarmatrixanimalmodel) that $G$ corresponds to the product of the numerator relationship matrix $A$ and the genetic additive variance $\sigma_u^2$. But the mixed model equations depend on the inverse $G^{-1}$ of $G$ which means
\begin{equation} G^{-1} = (A * \sigma_u^2)^{-1} = A^{-1} * {1\over \sigma_u^2} (#eq:ginvformme) \end{equation}
From \@ref(eq:ginvformme) we can see that in order to be able to set up the mixed model equations for an animal model, we need the inverse numerator relationship matrix $A^{-1}$. Because in practical routine predictions of breeding values, we have something in the order of 10 million animals that we predict breeding values for. Hence the matrix $A$ has 10 million rows and 10 million columns. A matrix of that size cannot be inverted explicitly with commonly known methods such as described under https://en.wikipedia.org/wiki/Invertible_matrix. This would have been the end of the BLUP animal model, if not Henderson published in [r mrmt$add('Henderson1976')] a set of simple rules to directly construct the matrix $A^{-1}$ without setting up the numerator relationship matrix.
Structure of $A^{-1}$
When looking at a concrete example of an inverse of a numerator relationship matrix as shown below, we can discover some important facts. Let us assume the following pedigree.
### # first six animals from Goetz p. 83 suppressPackageStartupMessages( library(pedigreemm) ) n_nr_ani_ped <- 5 n_nr_parent <- 3 tbl_ped_simple <- dplyr::data_frame(Calf = c(1:n_nr_ani_ped), Sire = c(NA, NA, NA, 1, 3), Dam = c(NA, NA, NA, 2, 2)) ped_simple <- pedigree(sire = tbl_ped_simple$Sire, dam = tbl_ped_simple$Dam, label = as.character(1:n_nr_ani_ped)) matu_simple <- as.matrix(getA(ped = ped_simple)) matAinv_simple <- as.matrix(getAInv(ped = ped_simple))
knitr::kable(tbl_ped_simple, format = 'latex', booktabs = TRUE, longtable = TRUE, caption = "Pedigree Used To Compute Inverse Numerator Relationship Matrix")
The numerator relationship matrix $A$ for the pedigree shown in Table \@ref(tab:simplepedexamplesetup) is shown in \@ref(eq:numeratorrelationshipmatrix).
cat("\\begin{equation}\n") # cat("A = \\left[") # cat(paste(rmddochelper::sConvertMatrixToLaTexArray(pmatAMatrix = matu_simple, pnDigits = 4), sep = "\n"), "\n") # cat("\\right]\n") cat(paste(rmdhelp::bmatrix(pmat = matu_simple, ps_name = 'A'), sep = '\n'), '\n') cat("\\label{eq:numeratorrelationshipmatrix}\n") cat("\\end{equation}\n")
The matrix $A^{-1}$ shown below corresponds to the inverse of the matrix $A$ from \@ref(eq:numeratorrelationshipmatrix).
cat("$$\n") # cat("A^{-1} = \\left[") # cat(paste(rmddochelper::sConvertMatrixToLaTexArray(pmatAMatrix = matAinv_simple, pnDigits = 4), sep = "\n"), "\n") # cat("\\right]") cat(paste(rmdhelp::bmatrix(pmat = matAinv_simple, ps_name = 'A^{-1}'), sep = '\n'), '\n') cat("$$")
From the above shown inverse $A^{-1}$, we recognize that $A^{-1}$ has a simpler structure than $A$ itself. The reason for this is that non-zero elements occur only at matrix elements of $A^{-1}$ corresponding to parents and progeny or to mates. Furthermore off-diagonal elements corresponding to mates are positive and elements corresponding to parents and progeny are negative. These observations were used by Henderson in [@Henderson1976] to come up with the rules described below.
Henderson's Rule To Set Up $A^{-1}$
We denote the row or column index corresponding to an animal of interest as $i$ and the row or column index belonging to the animals father as $s$ and the row or column index corresponding to animal $i$'s mother as $d$. The rules differentiate the following three cases
- both parents of animal $i$ are known
- only one parent of animal $i$ is known
- both parents of animal $i$ are unknown
Both Parents Known
- add $2$ to the diagonal-element $(i,i)$
- add $-1$ to off-diagonal elements $(s,i)$, $(i,s)$, $(d,i)$ and $(i,d)$
- add ${1\over 2}$ to elements $(s,s)$, $(d,d)$, $(s,d)$, $(d,s)$
Only One Parent Known
We assume that sire $s$ is known
- add ${4\over 3}$ to diagonal-element $(i,i)$
- add $-{2\over 3}$ to off-diagonal elements $(s,i)$, $(i,s)$
- add ${1\over 3}$ to element $(s,s)$
Both Parents Unknown
- add $1$ to diagonal-element $(i,i)$
The application of Henderson's rules to construct $A^{-1}$ directly will be a problem in one of the coming exercises. When applying the above described rules, it has to be noted well that this simple version of the rules are only valid for a pedigree without inbreeding. We will see in section \@ref(#derivationofhendersonsrules) how to extend the simple rules such that they can be used for a general pedigree with inbreeding.
Derivation of Henderson's Rules {#derivationofhendersonsrules}
Henderson's rules can be related to the so-called LDL-decomposition of the numerator relationship matrix $A$. The result of this decomposition consists of the equivalence between matrix $A$ and the product of three matrices $L$, $D$ and $L^T$. The matrix $L$ is a lower triangular matrix and the matrix $D$ is a diagonal matrix. The reason for why we are doing this decomposition of $A$ is that the matrices $L$ and $D$ can much easier be inverted than the matrix $A$. The LDL-decomposition is a general procedure in linear algebra and it can be applied to any symmetric and positive-definite matrix not only to numerator relationship matrices. But when the LDL-decomposition is applied to a numerator relationship matrix, the matrices $L$ and $D$ do also have a special genetic meaning. This meaning is revealed in the following derivation.
Decomposition of True Breeding Value and its Variance {#decompositionoftruebreedingvalueanditsvariance}
The true breeding value ($u_i$) of animal $i$ can be decomposed into the true breeding values of animal $i$'s parents $s$ and $d$ and the mendelian sampling term $m_i$
\begin{equation} u_i = {1\over 2} u_s + {1\over 2} u_d + m_i (#eq:decomposeanimaltruebreedingvalue) \end{equation}
Applying the decomposition shown in \@ref(eq:decomposeanimaltruebreedingvalue) for all animals in the pedigree and combining the decompositions into a matrix-vector notation, we get
\begin{equation} u = P \cdot u + m (#eq:decomposetruebreedingvaluematvec) \end{equation}
\begin{tabular}{lll} where & & \ & $u$ & vector of true breeding values for all animals \ & $P$ & matrix linking animals to their parents \ & $m$ & vector of mendelian sampling terms \end{tabular}
Equation \@ref(eq:decomposetruebreedingvaluematvec) can be interpreted as regression of the true breeding values onto parental breeding values. In such a regression the random term $m$ is like a residual term. The genetic meaning of $m$ corresponds to the deviation of $u_i$ from the full-sib average of the true breeding values within the nuclear family with parents $s$ and $d$. The term $m$ is called mendelian sampling term. The source of $m$ is the fact that during meiosis, parental alleles are randomly assigned to each progeny. Bulmer [r mrmt$add('Bulmer1971')] has shown that $m_i$ are independent from true breeding values of parents $s$ and $d$. Using this result, we take the variance on both sides of \@ref(eq:decomposeanimaltruebreedingvalue) leading to
\begin{align}
var(u_i) &= var({1\over 2} u_s + {1\over 2} u_d + m_i) \notag \
&= {1\over 4} var(u_s) + {1\over 4} var(u_d) + {1\over 2} cov(u_s, u_d) + var(m_i)
(#eq:vartruebreedingvalue)
\end{align}
From \@ref(eq:genvarcovarmatrixanimalmodel) together with the definition of the numerator relationship matrix $A$, we know that
\begin{align}
var(u_i) &= (1 + F_i) \sigma_u^2 \notag \
var(u_s) &= (1 + F_s) \sigma_u^2 \notag \
var(u_d) &= (1 + F_d) \sigma_u^2 \notag \
cov(u_s, u_d) &= (A)_{sd} \sigma_u^2 = 2F_i \sigma_u^2
(#eq:usedefnumeratorrelationshipmatrix)
\end{align}
Variance of Mendelian Sampling Terms {#variancemendeliansamplingterm}
Inserting the relations from \@ref(eq:usedefnumeratorrelationshipmatrix) into \@ref(eq:vartruebreedingvalue) and solving for $var(m_i)$ gives the following result
\begin{align}
var(m_i) &= var(u_i) - {1\over 4} var(u_s) - {1\over 4} var(u_d) - {1\over 2} cov(u_s, u_d) \notag \
&= (1 + F_i) \sigma_u^2 - {1\over 4}(1 + F_s) \sigma_u^2 - {1\over 4} (1 + F_d) \sigma_u^2 - {1\over 2} * 2 * F_i * \sigma_u^2 \notag \
&= \left({1\over 2} - {1\over 4}(F_s + F_d)\right) \sigma_u^2
(#eq:mendeliansamplingvariance)
\end{align}
In case where only parent $s$ of animal $i$ is known the terms in \@ref(eq:decomposeanimaltruebreedingvalue) and \@ref(eq:mendeliansamplingvariance) change to
\begin{equation} u_i = {1\over 2} u_s + {1\over 2} u_d + m_i (#eq:decomposeanimaltruebreedingvalue) \end{equation}
\begin{align}
u_i &= {1\over 2} u_s + m_i \notag \
var(m_i) &= \left(1 - {1\over 4}(1+ F_s)\right) \sigma_u^2 \notag \
&= \left({3\over 4} - {1\over 4}F_s\right) \sigma_u^2
(#eq:mendeliansamplingvarianceoneparentknown)
\end{align}
When both parents are unknown, we get
\begin{align}
u_i &= m_i \notag \
var(m_i) &= \sigma_u^2
(#eq:mendeliansamplingvariancebothparentunknown)
\end{align}
Decomposition of $A$ {#decompositionofa}
The true breeding values $u_s$ of sire $s$ and $u_d$ of dam $d$ can be decomposed in a similar way as shown for the true breeding value $u_i$ in \@ref(eq:decomposeanimaltruebreedingvalue).
\begin{align} u_s &= {1\over 2} u_{ss} + {1\over 2} u_{ds} + m_s \notag \ u_d &= {1\over 2} u_{sd} + {1\over 2} u_{dd} + m_d (#eq:decomposeparenttruebreedingvalue) \end{align}
\begin{tabular}{lll}
where & & \
& $ss$ & sire of $s$ \
& $ds$ & dam of $s$ \
& $sd$ & sire of $d$ \
& $dd$ & dam of $d$
\end{tabular}
Using \@ref(eq:decomposeparenttruebreedingvalue) together with \@ref(eq:decomposeanimaltruebreedingvalue) leads to the following expression for $u_i$
\begin{align} u_i &= {1\over 2} u_s + {1\over 2} u_d + m_i \notag \ &= {1\over 2} ({1\over 2} u_{ss} + {1\over 2} u_{ds} + m_s)+ {1\over 2} ({1\over 2} u_{sd} + {1\over 2} u_{dd} + m_d) + m_i \notag \ &= {1\over 4} u_{ss} + {1\over 4} u_{ds} + {1\over 4} u_{sd} + {1\over 4} u_{dd} + {1\over 2} m_s + {1\over 2} m_d + m_i \notag (#eq:recursivedecomposeanimaltruebreedingvalue) \end{align}
This type of decomposition can also be done for the grand-parents of animal $i$ and further back until we get to true breeding values of animals with unknown parents. Animals of ancestor generations with unknown parents are called founder population. The process of decomposing true breeding values back through all generations of a pedigree is called recursive decomposition of animal $i$'s true breeding value. Because we know from \@ref(eq:mendeliansamplingvariancebothparentunknown) that the decomposition of true breeding values $u_k$ for an animal $k$ with unknown parents is $u_k = m_k$, the result of the recursive decomposition of $u_i$ is a sum of mendelian sampling terms linking the ancestors of $i$ back to the founder population.
Ordering all animals in a pedigree according to their age and writing the result of the recursive decomposition of all true breeding values in matrix-vector notation leads to
\begin{equation} u = L \cdot m (#eq:resultrecursivedecomposition) \end{equation}
The matrix $L$ is a lower triangular matrix. The row corresponding to animal $i$ is the average of the rows in $L$ corresponding to parents $s$ and $d$ of $i$. The matrix $L$ contains the path of the gene flow from the base population to all animals in the population. From equation \@ref(eq:resultrecursivedecomposition), we are computing the variance of all true breeding values which leads to
\begin{align} var(u) &= var(L \cdot m) \notag \ &= L \cdot var(m) \cdot L^T \notag \ &= L \cdot (D \sigma_u^2) \cdot L^T \notag \ &= L \cdot D \cdot L^T \sigma_u^2 \notag &= A \sigma_u^2 (#eq:vararecursivedecomposition) \end{align}
From \@ref(eq:vararecursivedecomposition), the LDL-decomposition of $A$ follows directly as $A = LDL^T$. The single components $m_i$ are independent of each other. This also means that $cov(m_i, m_j) = 0$ for $i \ne j$. Hence the matrix $D$ is a diagonal matrix. In section \@ref(#variancemendeliansamplingterm) we have seen that $var(m_i)$ always contain $\sigma_u^2$ as a factor. Hence it is reasonable to express $var(m)$ as a product of a diagonal matrix $D$ times the genetic additive variance $\sigma_u^2$. The diagonal elements of matrix $D$ are computed as shown in section \@ref(#variancemendeliansamplingterm).
The mixed model equations require the inverse numerator relationship matrix $A^{-1}$. Thanks to the LDL-decomposition of $A$, we can compute $A^{-1}$ as
\begin{equation} A^{-1} = (L \cdot D \cdot L^T)^{-1} = (L^T)^{-1} \cdot D^{-1} \cdot L^{-1} (#eq:ainvldldecomp) \end{equation}
The inverse $D^{-1}$ is easy to compute, because $D$ is a diagonal matrix. As a consequence of that $D^{-1}$ is also a diagonal matrix where the elements $(D^{-1}){ii}$ correspond to the inverse of the elements of the original matrix $D$, i.e. $(D^{-1}){ii} = 1/(D)_{ii}$. The matrix $L^{-1}$ is also a lower triangular matrix and can easily computed based on the two relations for the vector $m$. Based on \@ref(eq:decomposetruebreedingvaluematvec), we know
\begin{equation} m = u - Pu = (I-P)u (#eq:mendeliansamplingiminuspa) \end{equation}
and from \@ref(eq:resultrecursivedecomposition), we get
\begin{equation} m = L^{-1}u (#eq:mendeliansamplinglinversea) \end{equation}
Setting both expressions for $m$ in \@ref(eq:mendeliansamplingiminuspa) and \@ref(eq:mendeliansamplinglinversea) equal can be used to obtain $L^{-1}$
\begin{equation} m = (I-P)u = L^{-1}u (#eq:linverseequation) \end{equation}
Therefore,
\begin{equation} L^{-1} = I - P (#eq:linverseresult) \end{equation}
General Version of Henderson's Rule
Based on the decomposition of $A^{-1}$ shown in \@ref(eq:ainvldldecomp), the general version of Henderson's rule where inbreeding of animals can be accounted for can be summarized as
- Start with a matrix $A^{-1}$ where all elements are set to $0$.
- Let $d^i$ be the $i$-th diagonal element of $D^{-1}$ for animal $i$, assuming $i$ has parents $s$ and $d$.
- Then add the following contributions to $A^{-1}$
- $d^i$ to the element $(i,i)$
- $-d^i/2$ to the elements $(s,i)$, $(i,s)$, $(d,i)$, $(i,d)$
- $d^i/4$ to the elements $(s,s)$, $(s,d)$, $(d,s)$, $(d,d)$
The contributions to rows and columns corresponding to parents $s$ and $d$ are only added, if they are known. Because the elements $d^i$ are dependent on the inbreeding coefficients $F_s$ and $F_d$ of parents $s$ and $d$, we have to find an efficient way to compute inbreeding coefficients for all animals in a population.
Computing Inbreeding Coefficients {#computinginbreedingcoefficients}
Inbreeding coefficients can be computed using different methods. From all these methods, we are just showing the one method which was described in [r mrmt$add('Quaas1976')]. This method is based on the properties of a second decomposition of the numerator relationship matrix $A$ which is called the cholesky decomposition. The cholesky decomposition of a matrix $A$ corresponds to
\begin{equation} A = RR^T (#eq:choleskymata) \end{equation}
where the matrix $R$ is a lower triangular matrix. At this point we have to be clear about this. In practical routine evaluations, we will not explicitly compute this decomposition, because we do not want to construct $A$ explicitly. We are just using the properties of \@ref(eq:choleskymata) to find an efficient way to compute the diagonal elements $(A){ii}$ of $A$ without constructing the complete matrix $A$. The diagonal elements $(A){ii}$ are important, because they contain the inbreeding coefficients ($F_i$), as $(A){ii} = 1 + F_i$. Based on \@ref(eq:choleskymata), $(A){ii}$ can be computed from the components of $R$ as
\begin{equation} (A){ii} = \sum{j=1}^i (R)_{ij}^2 (#eq:DiagElemMatA) \end{equation}
This can be shown with a small $3\times 3$ matrix $A$
nAnzAni <- 3 matA <- rmdhelp::matGetMatElem(psBaseElement = "(A)", pnNrRow = nAnzAni, pnNrCol = nAnzAni) matR <- rmdhelp::matLowerTri(psBaseElement = "(R)", pnNrRow = nAnzAni, pnNrCol = nAnzAni) cat("$$") cat(paste(rmdhelp::bmatrix(pmat = matA), collapse = "\n")) cat(" \n") cat(" = ") cat(paste(rmdhelp::bmatrix(pmat = matR), collapse = "\n")) cat(" \n") cat(" * ") cat(paste(rmdhelp::bmatrix(pmat = t(matR)), collapse = "\n")) cat("\n") cat("$$\n")
For the above shown example, the diagonal elements $(A)_{ii}$ are computed as
\begin{align} (A){11} &= (R){11}^2 \notag \ (A){22} &= (R){21}^2 + (R){22}^2 \notag \ (A){33} &= (R){31}^2 + (R){32}^2 + (R)_{33}^2 \notag (#eq:MatADiag) \end{align}
This shows that diagonal elements $(A)_{ii}$ can be computed using just all the components of row $i$ in $R$ up to the diagonal. Next, we have to show how to compute the components of $R$.
Recursive Computation of $R$
We are setting the two decompositions of $A$ equal which leads to
\begin{equation} A = R * R^T = L * D * L^T (#eq:CholEqLdl) \end{equation}
Let us write the matrix $R$ as a product of two matrices $L$ and $S$ where $L$ corresponds to the matrix from the LDL-decomposition and insert that product into \@ref(eq:CholEqLdl)
\begin{equation} A = R * R^T = L * S * (L * S)^T = L * S * S^T * L^T= L * D * L^T (#eq:CholRLSEqLdl) \end{equation}
From \@ref(eq:CholRLSEqLdl), we conclude that $D = S \cdot S^T$ where $S$ is also a diagonal matrix with elements $(S){ii} = \sqrt{(D){ii}}$. For our small example we get
\begin{equation}
matL <- rmdhelp::matLowerTri(psBaseElement = "(L)", pnNrRow = nAnzAni, pnNrCol = nAnzAni, pvecDiag = 1) matS <- rmdhelp::matDiag(psBaseElement = "(S)", pnNrRow = nAnzAni, pnNrCol = nAnzAni) cat(paste(rmdhelp::bmatrix(pmat = matR), collapse = "\n")) cat(" = \n") cat(paste(rmdhelp::bmatrix(pmat = matL), collapse = "\n")) cat(" * \n") cat(paste(rmdhelp::bmatrix(pmat = matS), collapse = "\n")) cat("\n")
(#eq:SmallExRLS) \end{equation}
From \@ref(eq:SmallExRLS), we can see that the diagonal elements $(R){ii}$ are equal to $(S){ii}$. Therefore
\begin{equation} (R){ii} = (S){ii} = \sqrt{(D)_{ii}} \label{eq:ComputeRii} \end{equation}
Earlier, we have seen that diagonal elements $(D)_{ii}$ of $D$ correspond to
\begin{equation} (D){ii} = {1\over 2}\ - {1\over 4}\ (F_s + F_d) = 1 - 0.25((A){ss} + (A)_{dd}) \label{eq:ComputeDii} \end{equation}
and hence
\begin{equation} (R){ii} = \sqrt{1 - 0.25((A){ss} + (A)_{dd})} \label{eq:ComputeRiiResult} \end{equation}
The components $(A){ss}$ and $(A){dd}$ correspond to diagonal elements of $A$ of parents of $s$ and $d$ which were computed earlier.
The off-diagonal elements of $R$ are computed as
\begin{equation} (R){ij} = (L){ij} * (S)_{jj} \label{eq:ComputeRij} \end{equation}
One property of the matrix $L$ is that any element $(L){ij}$ is equal to the average of elements $(L){sj}$ and $(L)_{dj}$, if $s$ and $d$ are parents of animal $i$. Using this we get
\begin{align} (R){ij} &= (L){ij} * (S){jj} \nonumber\ &= {1\over 2}((L){sj} + (L){dj}) * (S){jj} \nonumber\ &= {1\over 2}((R){sj} + (R){dj}) (#eq:ComputeRijRec) \end{align}
With that we have shown how to compute all elements of $R$ recursively. This requires an ordering of all animals according to their age.
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.