This script is a protocol of the demo solution that was during the exercise hour on 2019-03-18. The purpose was to introduce all the students to the concepts of using mixed linear effects models to predict breeding values using a BLUP animal model.
Read input from csv files, start with observations
s_phe_path <- "https://charlotte-ngs.github.io/GELASMSS2019/ex/w05/data_ex04_phe.csv"
tbl_phe <- readr::read_csv(file = s_phe_path)## Parsed with column specification:
## cols(
## Animal = col_double(),
## Observation = col_double()
## )Read the pedigree
s_ped_path <- "https://charlotte-ngs.github.io/GELASMSS2019/ex/w05/data_ex04_ped.csv"
tbl_ped <- readr::read_csv(file = s_ped_path)## Parsed with column specification:
## cols(
## Animal = col_double(),
## Sire = col_double(),
## Dam = col_double()
## )Create the design matrices X and Z
n_obs <- nrow(tbl_phe)
X <- matrix(1, nrow=n_obs, ncol=1)
n_animal <- nrow(tbl_ped)
Z <- diag(nrow = n_animal)pedigree based on tbl_ped
?pedigreemm::pedigreeconstruct pedigree object
ped <- pedigreemm::pedigree(sire = tbl_ped$Sire, dam = tbl_ped$Dam, label = as.character(tbl_ped$Animal))inverse of A
Ainv <- as.matrix(pedigreemm::getAInv(ped = ped))
Ainv## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
## 1 2.0 0.0 0.5 0 0.5 0.0 -1.0 0.0 -1.0 0 0 0 0 0 0 0 0 0 0
## 2 0.0 2.0 0.5 0 0.5 -1.0 0.0 -1.0 0.0 0 0 0 0 0 0 0 0 0 0
## 3 0.5 0.5 2.0 0 0.0 -1.0 -1.0 0.0 0.0 0 0 0 0 0 0 0 0 0 0
## 4 0.0 0.0 0.0 1 0.0 0.0 0.0 0.0 0.0 0 0 0 0 0 0 0 0 0 0
## 5 0.5 0.5 0.0 0 2.0 0.0 0.0 -1.0 -1.0 0 0 0 0 0 0 0 0 0 0
## 6 0.0 -1.0 -1.0 0 0.0 5.0 0.0 1.5 1.5 0 0 -1 0 0 -1 -1 -1 -1 0
## 7 -1.0 0.0 -1.0 0 0.0 0.0 5.0 2.5 0.5 -1 -1 0 -1 -1 0 0 0 0 -1
## 8 0.0 -1.0 0.0 0 -1.0 1.5 2.5 6.0 0.0 -1 -1 0 -1 0 -1 0 -1 -1 -1
## 9 -1.0 0.0 0.0 0 -1.0 1.5 0.5 0.0 4.0 0 0 -1 0 -1 0 -1 0 0 0
## 10 0.0 0.0 0.0 0 0.0 0.0 -1.0 -1.0 0.0 2 0 0 0 0 0 0 0 0 0
## 11 0.0 0.0 0.0 0 0.0 0.0 -1.0 -1.0 0.0 0 2 0 0 0 0 0 0 0 0
## 12 0.0 0.0 0.0 0 0.0 -1.0 0.0 0.0 -1.0 0 0 2 0 0 0 0 0 0 0
## 13 0.0 0.0 0.0 0 0.0 0.0 -1.0 -1.0 0.0 0 0 0 2 0 0 0 0 0 0
## 14 0.0 0.0 0.0 0 0.0 0.0 -1.0 0.0 -1.0 0 0 0 0 2 0 0 0 0 0
## 15 0.0 0.0 0.0 0 0.0 -1.0 0.0 -1.0 0.0 0 0 0 0 0 2 0 0 0 0
## 16 0.0 0.0 0.0 0 0.0 -1.0 0.0 0.0 -1.0 0 0 0 0 0 0 2 0 0 0
## 17 0.0 0.0 0.0 0 0.0 -1.0 0.0 -1.0 0.0 0 0 0 0 0 0 0 2 0 0
## 18 0.0 0.0 0.0 0 0.0 -1.0 0.0 -1.0 0.0 0 0 0 0 0 0 0 0 2 0
## 19 0.0 0.0 0.0 0 0.0 0.0 -1.0 -1.0 0.0 0 0 0 0 0 0 0 0 0 2
## 20 0.0 0.0 0.0 0 0.0 -1.0 0.0 0.0 -1.0 0 0 0 0 0 0 0 0 0 0
## 21 0.0 0.0 0.0 0 0.0 0.0 -1.0 -1.0 0.0 0 0 0 0 0 0 0 0 0 0
## 20 21
## 1 0 0
## 2 0 0
## 3 0 0
## 4 0 0
## 5 0 0
## 6 -1 0
## 7 0 -1
## 8 0 -1
## 9 -1 0
## 10 0 0
## 11 0 0
## 12 0 0
## 13 0 0
## 14 0 0
## 15 0 0
## 16 0 0
## 17 0 0
## 18 0 0
## 19 0 0
## 20 2 0
## 21 0 2lambda from Problem
genvar <- 25
resvar <- 75
lambda <- resvar/genvarvector of observations
y <- tbl_phe$Observation
y## [1] 100.430 103.396 114.458 100.068 104.144 117.524 97.744 111.926
## [9] 103.486 97.914 104.651 115.714 86.900 101.097 102.795 112.182
## [17] 109.295 105.271 91.744 101.132 107.385Matrix M corresponding to the coefficient matrix of the MME
xtx <- crossprod(X)
xtx## [,1]
## [1,] 21xtz <- crossprod(X,Z)
xtz## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 1 1 1 1 1 1 1 1 1 1 1 1 1
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
## [1,] 1 1 1 1 1 1 1 1ztx <- t(xtz)
ztx## [,1]
## [1,] 1
## [2,] 1
## [3,] 1
## [4,] 1
## [5,] 1
## [6,] 1
## [7,] 1
## [8,] 1
## [9,] 1
## [10,] 1
## [11,] 1
## [12,] 1
## [13,] 1
## [14,] 1
## [15,] 1
## [16,] 1
## [17,] 1
## [18,] 1
## [19,] 1
## [20,] 1
## [21,] 1ztzlainv <- crossprod(Z) + lambda * Ainv
ztzlainv## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
## 1 7.0 0.0 1.5 0 1.5 0.0 -3.0 0.0 -3.0 0 0 0 0 0 0 0 0 0 0
## 2 0.0 7.0 1.5 0 1.5 -3.0 0.0 -3.0 0.0 0 0 0 0 0 0 0 0 0 0
## 3 1.5 1.5 7.0 0 0.0 -3.0 -3.0 0.0 0.0 0 0 0 0 0 0 0 0 0 0
## 4 0.0 0.0 0.0 4 0.0 0.0 0.0 0.0 0.0 0 0 0 0 0 0 0 0 0 0
## 5 1.5 1.5 0.0 0 7.0 0.0 0.0 -3.0 -3.0 0 0 0 0 0 0 0 0 0 0
## 6 0.0 -3.0 -3.0 0 0.0 16.0 0.0 4.5 4.5 0 0 -3 0 0 -3 -3 -3 -3 0
## 7 -3.0 0.0 -3.0 0 0.0 0.0 16.0 7.5 1.5 -3 -3 0 -3 -3 0 0 0 0 -3
## 8 0.0 -3.0 0.0 0 -3.0 4.5 7.5 19.0 0.0 -3 -3 0 -3 0 -3 0 -3 -3 -3
## 9 -3.0 0.0 0.0 0 -3.0 4.5 1.5 0.0 13.0 0 0 -3 0 -3 0 -3 0 0 0
## 10 0.0 0.0 0.0 0 0.0 0.0 -3.0 -3.0 0.0 7 0 0 0 0 0 0 0 0 0
## 11 0.0 0.0 0.0 0 0.0 0.0 -3.0 -3.0 0.0 0 7 0 0 0 0 0 0 0 0
## 12 0.0 0.0 0.0 0 0.0 -3.0 0.0 0.0 -3.0 0 0 7 0 0 0 0 0 0 0
## 13 0.0 0.0 0.0 0 0.0 0.0 -3.0 -3.0 0.0 0 0 0 7 0 0 0 0 0 0
## 14 0.0 0.0 0.0 0 0.0 0.0 -3.0 0.0 -3.0 0 0 0 0 7 0 0 0 0 0
## 15 0.0 0.0 0.0 0 0.0 -3.0 0.0 -3.0 0.0 0 0 0 0 0 7 0 0 0 0
## 16 0.0 0.0 0.0 0 0.0 -3.0 0.0 0.0 -3.0 0 0 0 0 0 0 7 0 0 0
## 17 0.0 0.0 0.0 0 0.0 -3.0 0.0 -3.0 0.0 0 0 0 0 0 0 0 7 0 0
## 18 0.0 0.0 0.0 0 0.0 -3.0 0.0 -3.0 0.0 0 0 0 0 0 0 0 0 7 0
## 19 0.0 0.0 0.0 0 0.0 0.0 -3.0 -3.0 0.0 0 0 0 0 0 0 0 0 0 7
## 20 0.0 0.0 0.0 0 0.0 -3.0 0.0 0.0 -3.0 0 0 0 0 0 0 0 0 0 0
## 21 0.0 0.0 0.0 0 0.0 0.0 -3.0 -3.0 0.0 0 0 0 0 0 0 0 0 0 0
## 20 21
## 1 0 0
## 2 0 0
## 3 0 0
## 4 0 0
## 5 0 0
## 6 -3 0
## 7 0 -3
## 8 0 -3
## 9 -3 0
## 10 0 0
## 11 0 0
## 12 0 0
## 13 0 0
## 14 0 0
## 15 0 0
## 16 0 0
## 17 0 0
## 18 0 0
## 19 0 0
## 20 7 0
## 21 0 7Putting together all the components
M <- rbind(cbind(xtx,xtz), cbind(ztx, ztzlainv))
M## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
## 21 1.0 1.0 1.0 1 1.0 1.0 1.0 1.0 1.0 1 1 1 1 1 1 1 1 1
## 1 1 7.0 0.0 1.5 0 1.5 0.0 -3.0 0.0 -3.0 0 0 0 0 0 0 0 0 0
## 2 1 0.0 7.0 1.5 0 1.5 -3.0 0.0 -3.0 0.0 0 0 0 0 0 0 0 0 0
## 3 1 1.5 1.5 7.0 0 0.0 -3.0 -3.0 0.0 0.0 0 0 0 0 0 0 0 0 0
## 4 1 0.0 0.0 0.0 4 0.0 0.0 0.0 0.0 0.0 0 0 0 0 0 0 0 0 0
## 5 1 1.5 1.5 0.0 0 7.0 0.0 0.0 -3.0 -3.0 0 0 0 0 0 0 0 0 0
## 6 1 0.0 -3.0 -3.0 0 0.0 16.0 0.0 4.5 4.5 0 0 -3 0 0 -3 -3 -3 -3
## 7 1 -3.0 0.0 -3.0 0 0.0 0.0 16.0 7.5 1.5 -3 -3 0 -3 -3 0 0 0 0
## 8 1 0.0 -3.0 0.0 0 -3.0 4.5 7.5 19.0 0.0 -3 -3 0 -3 0 -3 0 -3 -3
## 9 1 -3.0 0.0 0.0 0 -3.0 4.5 1.5 0.0 13.0 0 0 -3 0 -3 0 -3 0 0
## 10 1 0.0 0.0 0.0 0 0.0 0.0 -3.0 -3.0 0.0 7 0 0 0 0 0 0 0 0
## 11 1 0.0 0.0 0.0 0 0.0 0.0 -3.0 -3.0 0.0 0 7 0 0 0 0 0 0 0
## 12 1 0.0 0.0 0.0 0 0.0 -3.0 0.0 0.0 -3.0 0 0 7 0 0 0 0 0 0
## 13 1 0.0 0.0 0.0 0 0.0 0.0 -3.0 -3.0 0.0 0 0 0 7 0 0 0 0 0
## 14 1 0.0 0.0 0.0 0 0.0 0.0 -3.0 0.0 -3.0 0 0 0 0 7 0 0 0 0
## 15 1 0.0 0.0 0.0 0 0.0 -3.0 0.0 -3.0 0.0 0 0 0 0 0 7 0 0 0
## 16 1 0.0 0.0 0.0 0 0.0 -3.0 0.0 0.0 -3.0 0 0 0 0 0 0 7 0 0
## 17 1 0.0 0.0 0.0 0 0.0 -3.0 0.0 -3.0 0.0 0 0 0 0 0 0 0 7 0
## 18 1 0.0 0.0 0.0 0 0.0 -3.0 0.0 -3.0 0.0 0 0 0 0 0 0 0 0 7
## 19 1 0.0 0.0 0.0 0 0.0 0.0 -3.0 -3.0 0.0 0 0 0 0 0 0 0 0 0
## 20 1 0.0 0.0 0.0 0 0.0 -3.0 0.0 0.0 -3.0 0 0 0 0 0 0 0 0 0
## 21 1 0.0 0.0 0.0 0 0.0 0.0 -3.0 -3.0 0.0 0 0 0 0 0 0 0 0 0
## 19 20 21
## 1 1 1
## 1 0 0 0
## 2 0 0 0
## 3 0 0 0
## 4 0 0 0
## 5 0 0 0
## 6 0 -3 0
## 7 -3 0 -3
## 8 -3 0 -3
## 9 0 -3 0
## 10 0 0 0
## 11 0 0 0
## 12 0 0 0
## 13 0 0 0
## 14 0 0 0
## 15 0 0 0
## 16 0 0 0
## 17 0 0 0
## 18 0 0 0
## 19 7 0 0
## 20 0 7 0
## 21 0 0 7defining the right-hand side
r <- rbind(crossprod(X,y), crossprod(Z,y))
r## [,1]
## [1,] 2189.256
## [2,] 100.430
## [3,] 103.396
## [4,] 114.458
## [5,] 100.068
## [6,] 104.144
## [7,] 117.524
## [8,] 97.744
## [9,] 111.926
## [10,] 103.486
## [11,] 97.914
## [12,] 104.651
## [13,] 115.714
## [14,] 86.900
## [15,] 101.097
## [16,] 102.795
## [17,] 112.182
## [18,] 109.295
## [19,] 105.271
## [20,] 91.744
## [21,] 101.132
## [22,] 107.385Computing the solution
l <- solve(M,r)
l## [,1]
## 104.19659668
## 1 -2.38439070
## 2 1.05346529
## 3 2.33905342
## 4 -1.03214917
## 5 0.02402116
## 6 3.97976622
## 7 -2.60790538
## 8 -0.07327769
## 9 -0.51860340
## 10 -2.04659227
## 11 -1.08416370
## 12 3.12869882
## 13 -3.62002084
## 14 -1.78273186
## 15 1.47398127
## 16 2.62412740
## 17 2.40255270
## 18 1.82769556
## 19 -2.92802084
## 20 1.04555597
## 21 -0.69359227