tbl_si_dat <- tibble::tibble(Animal = c(1:5),
Sire = c(NA, NA, NA, 1, 4),
Dam = c(NA, NA, NA, 2, 3),
Sex = c("M", "F", "F", "M", "F"),
WWG = c(4.5, 2.9, 3.9, 3.5, 5.0),
PWG = c(7, 4.2, 5.9, 6.1, 6.9))
tbl_si_datVector \(y\)
(vec_y <- c(tbl_si_dat$WWG, tbl_si_dat$PWG))
[1] 4.5 2.9 3.9 3.5 5.0 7.0 4.2 5.9 6.1 6.9The genetic variance-covariance matrix \(G_0\) contains variances and covariances between the different traits.
(mat_G0 <- matrix(c(20, 18,
18, 40), nrow = 2, byrow = TRUE))
[,1] [,2]
[1,] 20 18
[2,] 18 40This means that the genetic variance of trait 1 (‘WWG’) is 20 and the genetic variance of trait 2 (‘PWG’) is 40 and the genetic covariance between the two traits is 18.
Similarly for the errors, we the matrix \(R_0\)
(mat_R0 <- matrix(c(40, 11,
11, 30), nrow = 2, byrow = TRUE))
[,1] [,2]
[1,] 40 11
[2,] 11 30The variance covariance matrix of all the error terms (for all animals with observations and all traits),
\[R = R_0 \otimes I_n\]
(mat_R <- mat_R0 %x% mat_In)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 40 0 0 0 0 11 0 0 0 0
[2,] 0 40 0 0 0 0 11 0 0 0
[3,] 0 0 40 0 0 0 0 11 0 0
[4,] 0 0 0 40 0 0 0 0 11 0
[5,] 0 0 0 0 40 0 0 0 0 11
[6,] 11 0 0 0 0 30 0 0 0 0
[7,] 0 11 0 0 0 0 30 0 0 0
[8,] 0 0 11 0 0 0 0 30 0 0
[9,] 0 0 0 11 0 0 0 0 30 0
[10,] 0 0 0 0 11 0 0 0 0 30The variance-covariance matrix of all the breeding values
\[G = G_0 \otimes A\]
where \(A\) is the numerator relationship matrix
The numerator relationship matrix \(A\)
(mat_A <- as.matrix(pedigreemm::getA(ped = ped)))The matrix \(G\) is
(mat_G <- mat_G0 %x% mat_A)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 20.0 0.0 0 10 5.0 18.0 0.0 0 9 4.5
[2,] 0.0 20.0 0 10 5.0 0.0 18.0 0 9 4.5
[3,] 0.0 0.0 20 0 10.0 0.0 0.0 18 0 9.0
[4,] 10.0 10.0 0 20 10.0 9.0 9.0 0 18 9.0
[5,] 5.0 5.0 10 10 20.0 4.5 4.5 9 9 18.0
[6,] 18.0 0.0 0 9 4.5 40.0 0.0 0 20 10.0
[7,] 0.0 18.0 0 9 4.5 0.0 40.0 0 20 10.0
[8,] 0.0 0.0 18 0 9.0 0.0 0.0 40 0 20.0
[9,] 9.0 9.0 0 18 9.0 20.0 20.0 0 40 20.0
[10,] 4.5 4.5 9 9 18.0 10.0 10.0 20 20 40.0The matrix \(G^{-1} = G_0^{-1} \otimes A^{-1}\)
mat_G0inv <- solve(mat_G0)
mat_Ainv <- as.matrix(pedigreemm::getAInv(ped = ped))
(mat_Ginv <- mat_G0inv %x% mat_Ainv) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0.12605042 0.04201681 0.00000000 -0.08403361 0.00000000 -0.05672269 -0.01890756 0.00000000 0.03781513 0.00000000
[2,] 0.04201681 0.12605042 0.00000000 -0.08403361 0.00000000 -0.01890756 -0.05672269 0.00000000 0.03781513 0.00000000
[3,] 0.00000000 0.00000000 0.12605042 0.04201681 -0.08403361 0.00000000 0.00000000 -0.05672269 -0.01890756 0.03781513
[4,] -0.08403361 -0.08403361 0.04201681 0.21008403 -0.08403361 0.03781513 0.03781513 -0.01890756 -0.09453782 0.03781513
[5,] 0.00000000 0.00000000 -0.08403361 -0.08403361 0.16806723 0.00000000 0.00000000 0.03781513 0.03781513 -0.07563025
[6,] -0.05672269 -0.01890756 0.00000000 0.03781513 0.00000000 0.06302521 0.02100840 0.00000000 -0.04201681 0.00000000
[7,] -0.01890756 -0.05672269 0.00000000 0.03781513 0.00000000 0.02100840 0.06302521 0.00000000 -0.04201681 0.00000000
[8,] 0.00000000 0.00000000 -0.05672269 -0.01890756 0.03781513 0.00000000 0.00000000 0.06302521 0.02100840 -0.04201681
[9,] 0.03781513 0.03781513 -0.01890756 -0.09453782 0.03781513 -0.04201681 -0.04201681 0.02100840 0.10504202 -0.04201681
[10,] 0.00000000 0.00000000 0.03781513 0.03781513 -0.07563025 0.00000000 0.00000000 -0.04201681 -0.04201681 0.08403361The matrix \(G^{-1}\) is to be used in mixed model equations of Exercise 9.
round(mat_G %*% mat_Ginv, digits = 6)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 0 0 0 0 0 0 0 0 0
[2,] 0 1 0 0 0 0 0 0 0 0
[3,] 0 0 1 0 0 0 0 0 0 0
[4,] 0 0 0 1 0 0 0 0 0 0
[5,] 0 0 0 0 1 0 0 0 0 0
[6,] 0 0 0 0 0 1 0 0 0 0
[7,] 0 0 0 0 0 0 1 0 0 0
[8,] 0 0 0 0 0 0 0 1 0 0
[9,] 0 0 0 0 0 0 0 0 1 0
[10,] 0 0 0 0 0 0 0 0 0 1