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

Solutions

sol <- solve(coef_mat, rhs)
sol
         [,1]
   3.92136986
1  0.20091324
2 -0.33242009
3  0.13150685
4 -0.05369863
5  0.24684932

Problem 1: Decomposition

From the lecture, we know how the predicted breeding value of animal $4 can be decomposed into different components. Verify this decomposition based on the solutions obtained above.

Solution

Problem 2: Reliabilities

Compute the reliabilities of all the predicted breeding values.

Solution

LS0tCnRpdGxlOiBMaXZlc3RvY2sgQnJlZWRpbmcgYW5kIEdlbm9taWNzIC0gRXhlcmNpc2UgMTEKYXV0aG9yOiBQZXRlciB2b24gUm9ocgpkYXRlOiAiMjAyMS0xMi0wMyIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3Igc2V0dXAsIGluY2x1ZGU9RkFMU0V9CmtuaXRyOjpvcHRzX2NodW5rJHNldChlY2hvID0gVFJVRSkKYGBgCgojIERhdGEKCmBgYHtyfQpucl9hbmltYWwgPC0gNQp0YmxfZGVjb21wIDwtIHRpYmJsZTo6dGliYmxlKEFuaW1hbCA9IGMoMTpucl9hbmltYWwpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIFNpcmUgPSBjKE5BLCBOQSwgTkEsIDEsIDQpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgIERhbSA9IGMoTkEsIE5BLCBOQSwgMiwgMyksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgVHJhaXQgPSBjKDQuNSwgMi45LCAzLjksIDMuNSwgNS4wKSkKdGJsX2RlY29tcApgYGAKCgojIE1vZGVsCgokJHkgPSBYXG11ICsgWnUgKyBlJCQKCgojIENvbXBvbmVudHMKCmBgYHtyfQpYID0gbWF0cml4KDEsIG5yb3cgPSBucl9hbmltYWwsIG5jb2wgPSAxKQpYCmBgYAoKYGBge3J9ClogPSBkaWFnKDEsIG5yb3cgPSBucl9hbmltYWwpCloKYGBgCgoKIyBQZWRpZ3JlZQoKYGBge3J9CnBlZCA9IHBlZGlncmVlbW06OnBlZGlncmVlKHNpcmUgPSB0YmxfZGVjb21wJFNpcmUsIAogICAgICAgICAgICAgICAgICAgICAgICAgICBkYW0gPSB0YmxfZGVjb21wJERhbSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgbGFiZWwgPSBhcy5jaGFyYWN0ZXIoMTpucl9hbmltYWwpKQpBaW52IDwtIGFzLm1hdHJpeChwZWRpZ3JlZW1tOjpnZXRBSW52KHBlZCA9IHBlZCkpCkFpbnYKICAgICAgICAgICAgICAgICAgCmBgYAoKCiMgTU1FCgpgYGB7cn0KeHR4IDwtIGNyb3NzcHJvZChYKQp4dHgKYGBgCgpgYGB7cn0KeHR6IDwtIGNyb3NzcHJvZChYLCBaKQp4dHoKYGBgCgpgYGB7cn0KenR4IDwtIGNyb3NzcHJvZChaLCBYKQp6dHgKYGBgCgpgYGB7cn0KenR6IDwtIGNyb3NzcHJvZChaKQp6dHoKYGBgCgpgYGB7cn0Kc2lnbWFlMiA8LSA0MApzaWdtYXUyIDwtIDIwCmxhbWJkYSA8LSBzaWdtYWUyL3NpZ21hdTIKbGFtYmRhCmBgYAoKYGBge3J9Cnp0emFpbnZsYW1iZGEgPC0genR6ICsgQWludiAqIGxhbWJkYQp6dHphaW52bGFtYmRhCmBgYAoKYGBge3J9CmNvZWZfbWF0IDwtIHJiaW5kKGNiaW5kKHh0eCwgeHR6KSwgY2JpbmQoenR4LCB6dHphaW52bGFtYmRhKSkKY29lZl9tYXQKYGBgCgpgYGB7cn0KeSA8LSB0YmxfZGVjb21wJFRyYWl0Cnh0eSA8LSBjcm9zc3Byb2QoWCx5KQp6dHkgPC0gY3Jvc3Nwcm9kKFoseSkKcmhzIDwtIHJiaW5kKHh0eSwgenR5KQpyaHMgCgpgYGAKCiMgU29sdXRpb25zCgpgYGB7cn0Kc29sIDwtIHNvbHZlKGNvZWZfbWF0LCByaHMpCnNvbApgYGAKCgojIFByb2JsZW0gMTogRGVjb21wb3NpdGlvbgpGcm9tIHRoZSBsZWN0dXJlLCB3ZSBrbm93IGhvdyB0aGUgcHJlZGljdGVkIGJyZWVkaW5nIHZhbHVlIG9mIGFuaW1hbCAkNCBjYW4gYmUgZGVjb21wb3NlZCBpbnRvIGRpZmZlcmVudCBjb21wb25lbnRzLiBWZXJpZnkgdGhpcyBkZWNvbXBvc2l0aW9uIGJhc2VkIG9uIHRoZSBzb2x1dGlvbnMgb2J0YWluZWQgYWJvdmUuCgoKIyMgU29sdXRpb24KCgojIFByb2JsZW0gMjogUmVsaWFiaWxpdGllcwpDb21wdXRlIHRoZSByZWxpYWJpbGl0aWVzIG9mIGFsbCB0aGUgcHJlZGljdGVkIGJyZWVkaW5nIHZhbHVlcy4gCgoKIyMgU29sdXRpb24K