In charlotte-ngs/lbgfs2020: Livestock Breeding and Genomics (Fall 2020)

knitr::opts_chunk$set(echo = FALSE, results = "asis")
knitr::knit_hooks$set(hook_convert_odg = rmdhelp::hook_convert_odg)

Structure of $A^{-1}$

• Look at a simple example of $A$ and $A^{-1}$

### # first six animals from Goetz p. 83
suppressPackageStartupMessages( library(pedigreemm) )
n_nr_ani_ped <- 5
n_nr_parent <- 3
tbl_ped_simple <- tibble::tibble(Calf = c(1:n_nr_ani_ped),
                             Sire = c(NA, NA, NA, 1, 3),
                             Dam  = c(NA, NA, NA, 2, 2))
### # pedigreemm
ped_simple <- pedigree(sire = tbl_ped_simple$Sire, dam = tbl_ped_simple$Dam, label = as.character(1:n_nr_ani_ped))
matA_simple <- as.matrix(getA(ped = ped_simple))
matAinv_simple <- as.matrix(getAInv(ped = ped_simple))
### # LDL decomposition based on cholesky
matR <- t(chol(matA_simple))
### # matS = sqrt(matD)
matD <- diag(Dmat(ped = ped_simple), n_nr_ani_ped)
matS <- sqrt(matD)
matL <- matR %*% solve(matS)

knitr::kable(tbl_ped_simple,
             booktabs = TRUE,
             longtable = TRUE,
             caption = "Pedigree Used To Compute Inverse Numerator Relationship Matrix")

Numerator Relationship Matrix $A$

cat("\\begin{equation}\n")
cat("A = \\left[")
cat(paste(rmdhelp::sConvertMatrixToLaTexArray(pmatAMatrix = matA_simple, pnDigits = 4), sep = "\n"), "\n")
cat("\\right]\n")
cat("\\notag\n")
cat("\\end{equation}\n")

Inverse Numerator Relationship Matrix $A^{-1}$

cat("$$\n")
cat("A^{-1} = \\left[")
cat(paste(rmdhelp::sConvertMatrixToLaTexArray(pmatAMatrix = matAinv_simple, pnDigits = 4), sep = "\n"), "\n")
cat("\\right]")
cat("$$")

Conclusions

• $A^{-1}$ has simpler structure than $A$ itself
• Non-zero elements only at positions of parent-progeny and parent-mate positions
• Parent-mate positions are positive, parent-progeny are negative

Henderson's Rules

• Based on LDL-decomposition of $A$

$$A= L * D * L^T$$

\begin{tabular}{lll} where & $L$ & Lower triangular matrix \ & $D$ & Diagonal matrix \end{tabular}

• Why?
  • matrices $L$ and $D$ can be inverted directly, we 'll see how ...
  • construct $A^{-1} = (L^T)^{-1} * D^{-1} * L^{-1}$

Example

cat("$$\n")
cat("L = \\left[")
cat(paste(rmdhelp::sConvertMatrixToLaTexArray(pmatAMatrix = matL, pnDigits = 1), sep = "\n"), "\n")
cat("\\right]\n")
cat("$$\n")
cat("$$\n")
cat("D = \\left[")
cat(paste(rmdhelp::sConvertMatrixToLaTexArray(pmatAMatrix = matD, pnDigits = 1), sep = "\n"), "\n")
cat("\\right]\n")
cat("$$\n")

$\rightarrow$ Verify that $A = L * D * L^T$

Decomposition of True Breeding Value

• True breeding value ($u_i$) of animal $i$

$$ u_i = {1\over 2} u_s + {1\over 2} u_d + m_i $$

• Do that for all animals in pedigree

Decomposition for Example

\begin{align} u_1 &= m_1 \notag \ u_2 &= m_2 \notag \ u_3 &= m_3 \notag \ u_4 &= {1\over 2} u_1 + {1\over 2} u_2 + m_4 \notag \ u_5 &= {1\over 2} u_3 + {1\over 2} u_2 + m_5 \notag \end{align}

Matrix Vector Notation

• Define vectors $u$ and $m$ as
• Coefficients of $u_s$ and $u_d$ into matrix $P$

### # definition of vectors a and m
vec_u <- rmdhelp::vecGetVecElem(psBaseElement = "u", pnVecLen = n_nr_ani_ped, psResult = "latex")
vec_m <- rmdhelp::vecGetVecElem(psBaseElement = "m", pnVecLen = n_nr_ani_ped, psResult = "latex")
### # define matrix P
matP <- matrix(0, nrow = n_nr_ani_ped, ncol = n_nr_ani_ped)
matP[4,1] <- matP[4,2] <- matP[5,3] <- matP[5,2] <- 0.5
### # show both vectors
cat("$$\n")
cat("u = \\left[\n")
cat(paste(rmdhelp::sConvertVectorToLaTexArray(pvec_avector = vec_u), collapse = "\n"), "\n")
cat("\\right] \\text{, }\n")
cat("m = \\left[\n")
cat(paste(rmdhelp::sConvertVectorToLaTexArray(pvec_avector = vec_m), collapse = "\n"), "\n")
cat("\\right] \\text{, }\n")
cat("P = \\left[")
cat(paste(rmdhelp::sConvertMatrixToLaTexArray(pmatAMatrix = matP, pnDigits = 1), sep = "\n"), "\n")
cat("\\right]\n")
cat("$$\n")

• Result: Decomposition of true breeding values

$$u = P \cdot u + m$$

Decomposition of Variance

• Analogous decomposition of $var(u_i)$

\begin{align} var(u_i) &= var(1/2u_s + 1/2u_d + m_i) \notag \ &= var(1/2u_s) + var(1/2u_d) + {1\over 2} * cov(u_s, u_d) + var(m_i) \notag \ &= 1/4 var(u_s) + 1/4 var(u_d) + {1\over 2} * cov(u_s, u_d) + var(m_i) \notag \end{align}

• From the definition of $A$

\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 \notag \end{align}

Variance of Mendelian Sampling Terms

• What is $var(m_i)$?

• Solve equation for $var(u_i)$ for $var(m_i)$

\begin{align} var(m_i) &= var(u_i) - 1/4 var(u_s) - 1/4 var(u_d) - 2 * cov(u_s, u_d) \notag \end{align}

• Insert definitions from $A$

\begin{align} var(m_i) &= (1 + F_i) \sigma_u^2 - 1/4 (1 + F_s) \sigma_u^2 - 1/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 \notag \end{align}

• True, for both parents $s$ and $d$ of animal $i$ are known

Unknown Parents

• Only parent $s$ of animal $i$ is known

\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 \notag \end{align}

• Both parents are unknown

\begin{align} u_i &= m_i \notag \ var(m_i) &= \sigma_u^2 \notag \end{align}

Recursive Decomposition

• True breeding values of $s$ and $d$ can be decomposed into

\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 \notag \end{align}

\begin{tabular}{lll} where & $ss$ & sire of $s$ \ & $ds$ & dam of $s$ \ & $sd$ & sire of $d$ \ & $dd$ & dam of $d$
\end{tabular}

Example

• Add animal 6 with parents 4 and 5 to our example pedigree

n_nr_ani_ped <- 6
n_nr_parent <- 3
tbl_ped_ext <- tibble::tibble(Calf = c(1:n_nr_ani_ped),
                             Sire = c(NA, NA, NA, 1, 3, 4),
                             Dam  = c(NA, NA, NA, 2, 2, 5))
### # pedigreemm
ped_ext <- pedigree(sire = tbl_ped_ext$Sire, dam = tbl_ped_ext$Dam, label = as.character(1:n_nr_ani_ped))
matA_ext <- as.matrix(getA(ped = ped_ext))
matAinv_ext <- as.matrix(getAInv(ped = ped_ext))
### # LDL decomposition based on cholesky
matR_ext <- t(chol(matA_ext))
### # matS = sqrt(matD)
matD_ext <- diag(Dmat(ped = ped_ext), n_nr_ani_ped)
matS_ext <- sqrt(matD_ext)
matL_ext <- matR_ext %*% solve(matS_ext)
### # show table
knitr::kable(tbl_ped_ext, 
             booktabs  = TRUE,
             longtab   = TRUE)

First Step Of Decomposition

\begin{align} u_1 &= m_1 \notag \ u_2 &= m_2 \notag \ u_3 &= m_3 \notag \ u_4 &= {1\over 2} u_1 + {1\over 2} u_2 + m_4 \notag \ u_5 &= {1\over 2} u_3 + {1\over 2} u_2 + m_5 \notag \ u_6 &= {1\over 2} u_4 + {1\over 2} u_5 + m_6 \notag \end{align}

Decompose Parents

\begin{align} u_1 &= m_1 \notag \ u_2 &= m_2 \notag \ u_3 &= m_3 \notag \ u_4 &= {1\over 2} m_1 + {1\over 2} m_2 + m_4 \notag \ u_5 &= {1\over 2} m_3 + {1\over 2} m_2 + m_5 \notag \ u_6 &= {1\over 2} \left({1\over 2}(u_1 + u_2) + m_4\right) + {1\over 2} \left({1\over 2} (u_3 + u_2) + m_5 \right) + m_6 \notag \ &= {1\over 4} (u_1 + u_2) + {1\over 2} m_4 + {1\over 4} (u_3 + u_2) + {1\over 2} m_5 + m_6 \notag \end{align}

Decompose Grand Parents

• Only animal 6 has true breeding values for grand parents

\begin{align} u_6 &= {1\over 4} (u_1 + u_2) + {1\over 2} m_4 + {1\over 4} (u_3 + u_2) + {1\over 2} m_5 + m_6 \notag \ &= {1\over 4} m_1 + {1\over 4} m_2 + {1\over 4} m_3 + {1\over 4} m_2 + {1\over 2} m_4 + {1\over 2} m_5 + m_6 \notag \ &= {1\over 4} m_1 + {1\over 2} m_2 + {1\over 4} m_3 + {1\over 2} m_4 + {1\over 2} m_5 + m_6 \notag \end{align}

Summary

\begin{align} u_1 &= m_1 \notag \ u_2 &= m_2 \notag \ u_3 &= m_3 \notag \ u_4 &= {1\over 2} m_1 + {1\over 2} m_2 + m_4 \notag \ u_5 &= {1\over 2} m_3 + {1\over 2} m_2 + m_5 \notag \ u_6 &= {1\over 4} m_1 + {1\over 2} m_2 + {1\over 4} m_3 + {1\over 2} m_4 + {1\over 2} m_5 + m_6 \notag \end{align}

Matrix-Vector Notation

• Use vectors $u$ and $m$ again

### # definition of vectors a and m
vec_u_ext <- rmdhelp::vecGetVecElem(psBaseElement = "u", pnVecLen = n_nr_ani_ped, psResult = "latex")
vec_m_ext <- rmdhelp::vecGetVecElem(psBaseElement = "m", pnVecLen = n_nr_ani_ped, psResult = "latex")
### # show both vectors
cat("\\small \n $$\n")
cat("u = \\left[\n")
cat(paste(rmdhelp::sConvertVectorToLaTexArray(pvec_avector = vec_u_ext), collapse = "\n"), "\n")
cat("\\right] \\text{, }\n")
cat("m = \\left[\n")
cat(paste(rmdhelp::sConvertVectorToLaTexArray(pvec_avector = vec_m_ext), collapse = "\n"), "\n")
cat("\\right] \\text{, }\n")
cat("L = \\left[")
cat(paste(rmdhelp::sConvertMatrixToLaTexArray(pmatAMatrix = matL_ext, pnDigits = 2), sep = "\n"), "\n")
cat("\\right]\n")
cat("$$\n \\normalsize \n")

• Result of recursive decomposition of $u_i$

$$u = L \cdot m$$

Property of $L$

• Meaning of Element $(L)_{ij}$ of Matrix $L$:

#rmdhelp::use_odg_graphic(ps_path = 'odg/property-matL.odg')
knitr::include_graphics(path = "odg/property-matL.png")

Property of $L$ II

• Element $(L)_{ij}$ ($i>j$) is the proportion of $m_j$ in $u_i$
• Given: $i$ has parents $s$ and $d$
• $m_j$ can only come from $u_s$ and $u_d$, because $u_i = 1/2 u_s + 1/2 u_d + m_i$
• The proportion of $m_j$ in $u_i$ is half the proportion of $m_j$ in $u_s$ and half the proportion of $m_j$ in $u_d$

$$\rightarrow L_{ij} = {1\over 2}L_{sj} + {1\over 2}L_{dj}$$

Example

• $L_{41}$, $L_{62}$

#rmdhelp::use_odg_graphic(ps_path = 'odg/ex-matL.odg')
knitr::include_graphics(path = "odg/ex-matL.png")

Variance From Recursive Decomposition

\begin{align} var(u) &= var(L \cdot m) \notag \ &= L \cdot var(m) \cdot L^T \notag \end{align}

where $var(m)$ is the variance-covariance matrix of all components in vector $m$.

• covariances of components $m_i$, $cov(m_i, m_j) = 0$ for $i \ne j$
• $var(m_i)$ computed as shown before

Result

• variance-covariance matrix $var(m)$ can be written as $D * \sigma_u^2$ where $D$ is diagnoal

\begin{align} \rightarrow var(u) &= 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 \notag \end{align}

\begin{align} \rightarrow A &= L \cdot D\cdot L^T \notag \end{align}

Inverse of $A$ Based on $L$ and $D$

• Matrix $A$ was decomposed into $A = L \cdot D \cdot L^T$
• Get $A^{-1}$ as $A^{-1} = (L^T)^{-1} D^{-1} L^{-1}$
• $D^{-1}$ is diagonal again with elements

$$(D^{-1}){ii} = 1/(D){ii}$$

Inverse of $L$

• Compute $m$ based on the two decompositions of $u$

$$u = P \cdot u + m \quad \text{and} \quad u = L \cdot m$$

• Solve both for $m$ and set them equal

$$m = u - P \cdot u = (I-P)\cdot u \quad \text{and} \quad m = L^{-1} \cdot u$$

$$(I-P) \cdot u = L^{-1} \cdot u$$ and

$$L^{-1} = I-P$$

Example

### # 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))
### # pedigreemm
ped_simple <- pedigree(sire = tbl_ped_simple$Sire, dam = tbl_ped_simple$Dam, label = as.character(1:n_nr_ani_ped))
matA_simple <- as.matrix(getA(ped = ped_simple))
matAinv_simple <- as.matrix(getAInv(ped = ped_simple))
### # LDL decomposition based on cholesky
matR <- t(chol(matA_simple))
### # matS = sqrt(matD)
matD <- diag(Dmat(ped = ped_simple), n_nr_ani_ped)
matS <- sqrt(matD)
matL <- matR %*% solve(matS)
### # inverse of D
matDinv <- solve(matD)
### # matrix P
matP <- diag(1, n_nr_ani_ped) - solve(matL)
matLinv <- diag(1, n_nr_ani_ped) - matP
### # decomposition of AInv
matLinvtmatDinv <- crossprod(matLinv, matDinv)

knitr::kable(tbl_ped_simple,
             booktabs = TRUE,
             longtable = TRUE)

Matrix $D^{-1}$

• Because $D$ is diagonal

cat("$$\n")
cat("D = \\left[")
cat(paste(rmdhelp::sConvertMatrixToLaTexArray(pmatAMatrix = matD, pnDigits = 1), sep = "\n"), "\n")
cat("\\right]\n")
cat("$$\n")

• We get $D^{-1}$ as

cat("$$\n")
cat("D^{-1} = \\left[")
cat(paste(rmdhelp::sConvertMatrixToLaTexArray(pmatAMatrix = matDinv, pnDigits = 1), sep = "\n"), "\n")
cat("\\right]\n")
cat("$$\n")

Matrix $L^{-1}$

• Use $L^{-1} = I-P$
• Matrix $P$ from simple decomposition

cat("$$\n")
cat("P = \\left[")
cat(paste(rmdhelp::sConvertMatrixToLaTexArray(pmatAMatrix = matP, pnDigits = 1), sep = "\n"), "\n")
cat("\\right]\n")
cat("$$\n")

• Therefore

cat("$$\n")
cat("L^{-1} = I-P = \\left[")
cat(paste(rmdhelp::sConvertMatrixToLaTexArray(pmatAMatrix = matLinv, pnDigits = 1), sep = "\n"), "\n")
cat("\\right]\n")
cat("$$\n")

Decomposition of $A^{-1}$ I

cat("$$\n")
cat("A^{-1} = (L^{-1})^T \\cdot D^{-1} \\cdot L^{-1}")
cat("$$\n")

cat("$$\n")
cat("(L^{-1})^T \\cdot D^{-1}")
cat("$$\n")

cat("\\small \n $$\n")
cat("\\left[")
cat(paste(rmdhelp::sConvertMatrixToLaTexArray(pmatAMatrix = t(matLinv), pnDigits = 1), sep = "\n"), "\n")
cat("\\right] \\cdot \n")
cat("\\left[")
cat(paste(rmdhelp::sConvertMatrixToLaTexArray(pmatAMatrix = matDinv, pnDigits = 1), sep = "\n"), "\n")
cat("\\right]\n")
cat("$$\n")
cat("$$ = \n")
cat("\\left[")
cat(paste(rmdhelp::sConvertMatrixToLaTexArray(pmatAMatrix = matLinvtmatDinv, pnDigits = 1), sep = "\n"), "\n")
cat("\\right]\n")
cat("$$ \\normalsize \n")

Decomposition of $A^{-1}$ II

cat("$$\n")
cat("A^{-1} = (L^{-1})^T \\cdot D^{-1} \\cdot L^{-1}")
cat("$$\n")

cat("\\small \n $$\n")
cat("\\left[")
cat(paste(rmdhelp::sConvertMatrixToLaTexArray(pmatAMatrix = matLinvtmatDinv, pnDigits = 1), sep = "\n"), "\n")
cat("\\right] \\cdot \n")
cat("\\left[")
cat(paste(rmdhelp::sConvertMatrixToLaTexArray(pmatAMatrix = matLinv, pnDigits = 1), sep = "\n"), "\n")
cat("\\right]\n")
cat("$$\n")
cat("$$ = \n")
cat("\\left[")
cat(paste(rmdhelp::sConvertMatrixToLaTexArray(pmatAMatrix = crossprod(t(matLinvtmatDinv),matLinv), pnDigits = 1), sep = "\n"), "\n")
cat("\\right]\n")
cat("$$ \\normalsize \n")

Henderson's Rules

• 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

  • 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)$

• Valid without inbreeding

charlotte-ngs/lbgfs2020 documentation built on Dec. 20, 2020, 5:39 p.m.