decomposition of Predicted BLUP Breeding Values

Data

nr_animal <- 5
tbl_decomp <- tibble::tibble(Animal = c(1:nr_animal),
                             Sire = c(NA, NA, NA, 1, 4),
                             Dam = c(NA, NA, NA, 2, 3),
                             Trait = c(4.5, 2.9, 3.9, 3.5, 5.0))
tbl_decomp

Model

\[y = X\mu + Zu + e\]

Components

X = matrix(1, nrow = nr_animal, ncol = 1)
X

     [,1]
[1,]    1
[2,]    1
[3,]    1
[4,]    1
[5,]    1

Z = diag(1, nrow = nr_animal)
Z

     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    0    0
[2,]    0    1    0    0    0
[3,]    0    0    1    0    0
[4,]    0    0    0    1    0
[5,]    0    0    0    0    1

Pedigree

ped = pedigreemm::pedigree(sire = tbl_decomp$Sire, 
                           dam = tbl_decomp$Dam,
                           label = as.character(1:nr_animal))
Ainv <- as.matrix(pedigreemm::getAInv(ped = ped))
Ainv

     1    2    3    4  5
1  1.5  0.5  0.0 -1.0  0
2  0.5  1.5  0.0 -1.0  0
3  0.0  0.0  1.5  0.5 -1
4 -1.0 -1.0  0.5  2.5 -1
5  0.0  0.0 -1.0 -1.0  2

MME

xtx <- crossprod(X)
xtx

     [,1]
[1,]    5

xtz <- crossprod(X, Z)
xtz

     [,1] [,2] [,3] [,4] [,5]
[1,]    1    1    1    1    1

ztx <- crossprod(Z, X)
ztx

     [,1]
[1,]    1
[2,]    1
[3,]    1
[4,]    1
[5,]    1

ztz <- crossprod(Z)
ztz

     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    0    0
[2,]    0    1    0    0    0
[3,]    0    0    1    0    0
[4,]    0    0    0    1    0
[5,]    0    0    0    0    1

sigmae2 <- 40
sigmau2 <- 20
lambda <- sigmae2/sigmau2
lambda

[1] 2

ztzainvlambda <- ztz + Ainv * lambda
ztzainvlambda

   1  2  3  4  5
1  4  1  0 -2  0
2  1  4  0 -2  0
3  0  0  4  1 -2
4 -2 -2  1  6 -2
5  0  0 -2 -2  5

coef_mat <- rbind(cbind(xtx, xtz), cbind(ztx, ztzainvlambda))
coef_mat

     1  2  3  4  5
  5  1  1  1  1  1
1 1  4  1  0 -2  0
2 1  1  4  0 -2  0
3 1  0  0  4  1 -2
4 1 -2 -2  1  6 -2
5 1  0  0 -2 -2  5

y <- tbl_decomp$Trait
xty <- crossprod(X,y)
zty <- crossprod(Z,y)
rhs <- rbind(xty, zty)
rhs 

     [,1]
[1,] 19.8
[2,]  4.5
[3,]  2.9
[4,]  3.9
[5,]  3.5
[6,]  5.0

Solitions

sol <- solve(coef_mat, rhs)
sol

         [,1]
   3.92136986
1  0.20091324
2 -0.33242009
3  0.13150685
4 -0.05369863
5  0.24684932

Decomposition