# The Repeatability Model - Date : 2020-11-11
# An example of a single diary herd
# Linear Models for the Prediction of Animal Breeding Values, Raphael Mrode, 3rd Edition
# Julius Mugambe(PhD)
# Clear works pace
rm(list = ls())
# Set working directory
# Setting working directory
setwd("G:/R-MME")
getwd()
## [1] "G:/R-MME"
# install.packages('MASS', dependencies = TRUE)
require(MASS)
## Loading required package: MASS
#inputting variance components
sigma_2_a = 20
sigma_2_e = 28
sigma_2_pe = 12
sigma_2_y = 60
alpha_1 = sigma_2_e/sigma_2_a
alpha_1
## [1] 1.4
alpha_2 = sigma_2_e/sigma_2_pe
alpha_2
## [1] 2.333333
# Repeatability (r)
repeatability = (sigma_2_a + sigma_2_pe)/sigma_2_y
repeatability
## [1] 0.5333333
# Load data
data <- read.csv("data_dairy.txt", sep="")
data
## Cow Sire Dam Parity HYS Fatyield.kg.
## 1 4 1 2 1 1 201
## 2 4 1 2 2 3 280
## 3 5 3 2 1 1 150
## 4 5 3 2 2 4 200
## 5 6 1 5 1 2 160
## 6 6 1 5 2 3 190
## 7 7 3 4 1 1 180
## 8 7 3 4 2 3 250
## 9 8 1 7 1 2 285
## 10 8 1 7 2 4 300
# Making design matrices using the function desgn()
# make design matrix
desgn <- function(v) {
if (is.numeric(v)) {
va = v
mrow = length(va)
mcol = max(va)
}
if (is.character(v)) {
vf = factor(v)
# Save the index value of each level
va = as.numeric(vf)
mrow = length(va)
mcol = length(levels(vf))
}
# Initializing X-matrix
X = matrix(data = c(0), nrow = mrow, ncol = mcol)
for (i in 1:mrow) {
ic = va[i]
X[i, ic] = 1
}
return(X)
}
# X matrix related the records of HYS and parity
# design matrix for HYS
x1 = desgn(data$HYS)
x1
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 0 0 1 0
## [3,] 1 0 0 0
## [4,] 0 0 0 1
## [5,] 0 1 0 0
## [6,] 0 0 1 0
## [7,] 1 0 0 0
## [8,] 0 0 1 0
## [9,] 0 1 0 0
## [10,] 0 0 0 1
# design matrix for parity
x2 = desgn(data$Parity)
x2
## [,1] [,2]
## [1,] 1 0
## [2,] 0 1
## [3,] 1 0
## [4,] 0 1
## [5,] 1 0
## [6,] 0 1
## [7,] 1 0
## [8,] 0 1
## [9,] 1 0
## [10,] 0 1
# Overall X matrix
X = cbind(x1,x2)
X
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1 0 0 0 1 0
## [2,] 0 0 1 0 0 1
## [3,] 1 0 0 0 1 0
## [4,] 0 0 0 1 0 1
## [5,] 0 1 0 0 1 0
## [6,] 0 0 1 0 0 1
## [7,] 1 0 0 0 1 0
## [8,] 0 0 1 0 0 1
## [9,] 0 1 0 0 1 0
## [10,] 0 0 0 1 0 1
# Z matrix for only the animals with records
Z = desgn(data$Cow)
Z
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 0 0 0 1 0 0 0 0
## [2,] 0 0 0 1 0 0 0 0
## [3,] 0 0 0 0 1 0 0 0
## [4,] 0 0 0 0 1 0 0 0
## [5,] 0 0 0 0 0 1 0 0
## [6,] 0 0 0 0 0 1 0 0
## [7,] 0 0 0 0 0 0 1 0
## [8,] 0 0 0 0 0 0 1 0
## [9,] 0 0 0 0 0 0 0 1
## [10,] 0 0 0 0 0 0 0 1
# W ( In this example Z and W are the same or equal)
W = Z
W
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 0 0 0 1 0 0 0 0
## [2,] 0 0 0 1 0 0 0 0
## [3,] 0 0 0 0 1 0 0 0
## [4,] 0 0 0 0 1 0 0 0
## [5,] 0 0 0 0 0 1 0 0
## [6,] 0 0 0 0 0 1 0 0
## [7,] 0 0 0 0 0 0 1 0
## [8,] 0 0 0 0 0 0 1 0
## [9,] 0 0 0 0 0 0 0 1
## [10,] 0 0 0 0 0 0 0 1
# Diagonal term for the Incidence matrix for all the animals with records
I = diag((8))
I
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1 0 0 0 0 0 0 0
## [2,] 0 1 0 0 0 0 0 0
## [3,] 0 0 1 0 0 0 0 0
## [4,] 0 0 0 1 0 0 0 0
## [5,] 0 0 0 0 1 0 0 0
## [6,] 0 0 0 0 0 1 0 0
## [7,] 0 0 0 0 0 0 1 0
## [8,] 0 0 0 0 0 0 0 1
# y (phenotypes)
y = data$Fatyield.kg.
y
## [1] 201 280 150 200 160 190 180 250 285 300
# load pedigree data
ped <- read.csv("pedigree_dairy.txt", sep="")
# function to make inverse of numerator relationship matrix
ainv = function(ped) {
n = nrow(ped)
Ainv = matrix(c(0), nrow = n, ncol = n)
for (i in 1:n) {
animal = ped[i, 1]
sire = ped[i, 2]
dam = ped[i, 3]
if (sire == 0 & dam == 0) {
# both parents unknown
alpha = 1
Ainv[animal, animal] = alpha + Ainv[animal, animal]
} else if (sire != 0 & dam == 0) {
# sire known
alpha = 4/3
Ainv[animal, animal] = alpha + Ainv[animal, animal]
Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
} else if (sire == 0 & dam != 0) {
# dam known
alpha = 4/3
Ainv[animal, animal] = alpha + Ainv[animal, animal]
Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
} else {
# both parents known
alpha = 2
Ainv[animal, animal] = alpha + Ainv[animal, animal]
Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
Ainv[sire, dam] = alpha/4 + Ainv[sire, dam]
Ainv[dam, sire] = alpha/4 + Ainv[dam, sire]
Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
}
}
return(Ainv)
}
# Calculate the AINV from the pedigree using the ainv function
A <- ainv(ped)
A
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 2.5 0.5 0.0 -1.0 0.5 -1 0.5 -1
## [2,] 0.5 2.0 0.5 -1.0 -1.0 0 0.0 0
## [3,] 0.0 0.5 2.0 0.5 -1.0 0 -1.0 0
## [4,] -1.0 -1.0 0.5 2.5 0.0 0 -1.0 0
## [5,] 0.5 -1.0 -1.0 0.0 2.5 -1 0.0 0
## [6,] -1.0 0.0 0.0 0.0 -1.0 2 0.0 0
## [7,] 0.5 0.0 -1.0 -1.0 0.0 0 2.5 -1
## [8,] -1.0 0.0 0.0 0.0 0.0 0 -1.0 2
# LHS, RHS
# LHS construction
LHS = rbind(
cbind(t(X) %*% X, t(X) %*% Z, t(X) %*% W),
cbind(t(Z) %*% X, t(Z) %*% Z + A * alpha_1, t(Z) %*% W),
cbind(t(W) %*% X, t(W) %*% Z, t(W) %*% W + I * alpha_2))
# print
LHS
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 3 0 0 0 3 0 0.0 0.0 0.0 1.0 1.0 0.0 1.0
## [2,] 0 2 0 0 2 0 0.0 0.0 0.0 0.0 0.0 1.0 0.0
## [3,] 0 0 3 0 0 3 0.0 0.0 0.0 1.0 0.0 1.0 1.0
## [4,] 0 0 0 2 0 2 0.0 0.0 0.0 0.0 1.0 0.0 0.0
## [5,] 3 2 0 0 5 0 0.0 0.0 0.0 1.0 1.0 1.0 1.0
## [6,] 0 0 3 2 0 5 0.0 0.0 0.0 1.0 1.0 1.0 1.0
## [7,] 0 0 0 0 0 0 3.5 0.7 0.0 -1.4 0.7 -1.4 0.7
## [8,] 0 0 0 0 0 0 0.7 2.8 0.7 -1.4 -1.4 0.0 0.0
## [9,] 0 0 0 0 0 0 0.0 0.7 2.8 0.7 -1.4 0.0 -1.4
## [10,] 1 0 1 0 1 1 -1.4 -1.4 0.7 5.5 0.0 0.0 -1.4
## [11,] 1 0 0 1 1 1 0.7 -1.4 -1.4 0.0 5.5 -1.4 0.0
## [12,] 0 1 1 0 1 1 -1.4 0.0 0.0 0.0 -1.4 4.8 0.0
## [13,] 1 0 1 0 1 1 0.7 0.0 -1.4 -1.4 0.0 0.0 5.5
## [14,] 0 1 0 1 1 1 -1.4 0.0 0.0 0.0 0.0 0.0 -1.4
## [15,] 0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## [16,] 0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## [17,] 0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## [18,] 1 0 1 0 1 1 0.0 0.0 0.0 2.0 0.0 0.0 0.0
## [19,] 1 0 0 1 1 1 0.0 0.0 0.0 0.0 2.0 0.0 0.0
## [20,] 0 1 1 0 1 1 0.0 0.0 0.0 0.0 0.0 2.0 0.0
## [21,] 1 0 1 0 1 1 0.0 0.0 0.0 0.0 0.0 0.0 2.0
## [22,] 0 1 0 1 1 1 0.0 0.0 0.0 0.0 0.0 0.0 0.0
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
## [1,] 0.0 0.000000 0.000000 0.000000 1.000000 1.000000 0.000000 1.000000
## [2,] 1.0 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000
## [3,] 0.0 0.000000 0.000000 0.000000 1.000000 0.000000 1.000000 1.000000
## [4,] 1.0 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000 0.000000
## [5,] 1.0 0.000000 0.000000 0.000000 1.000000 1.000000 1.000000 1.000000
## [6,] 1.0 0.000000 0.000000 0.000000 1.000000 1.000000 1.000000 1.000000
## [7,] -1.4 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
## [8,] 0.0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
## [9,] 0.0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
## [10,] 0.0 0.000000 0.000000 0.000000 2.000000 0.000000 0.000000 0.000000
## [11,] 0.0 0.000000 0.000000 0.000000 0.000000 2.000000 0.000000 0.000000
## [12,] 0.0 0.000000 0.000000 0.000000 0.000000 0.000000 2.000000 0.000000
## [13,] -1.4 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.000000
## [14,] 4.8 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
## [15,] 0.0 2.333333 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
## [16,] 0.0 0.000000 2.333333 0.000000 0.000000 0.000000 0.000000 0.000000
## [17,] 0.0 0.000000 0.000000 2.333333 0.000000 0.000000 0.000000 0.000000
## [18,] 0.0 0.000000 0.000000 0.000000 4.333333 0.000000 0.000000 0.000000
## [19,] 0.0 0.000000 0.000000 0.000000 0.000000 4.333333 0.000000 0.000000
## [20,] 0.0 0.000000 0.000000 0.000000 0.000000 0.000000 4.333333 0.000000
## [21,] 0.0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 4.333333
## [22,] 2.0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
## [,22]
## [1,] 0.000000
## [2,] 1.000000
## [3,] 0.000000
## [4,] 1.000000
## [5,] 1.000000
## [6,] 1.000000
## [7,] 0.000000
## [8,] 0.000000
## [9,] 0.000000
## [10,] 0.000000
## [11,] 0.000000
## [12,] 0.000000
## [13,] 0.000000
## [14,] 2.000000
## [15,] 0.000000
## [16,] 0.000000
## [17,] 0.000000
## [18,] 0.000000
## [19,] 0.000000
## [20,] 0.000000
## [21,] 0.000000
## [22,] 4.333333
# RHS construction
RHS = rbind(t(X) %*% y, t(Z) %*% y, t(W) %*% y)
# print
RHS
## [,1]
## [1,] 531
## [2,] 445
## [3,] 720
## [4,] 500
## [5,] 976
## [6,] 1220
## [7,] 0
## [8,] 0
## [9,] 0
## [10,] 481
## [11,] 350
## [12,] 350
## [13,] 430
## [14,] 585
## [15,] 0
## [16,] 0
## [17,] 0
## [18,] 481
## [19,] 350
## [20,] 350
## [21,] 430
## [22,] 585
# Solutions
# generalized inverse of LHS
gi_LHS = ginv(LHS)
# print
gi_LHS
## [,1] [,2] [,3] [,4] [,5]
## [1,] 3.362958e-01 -2.302930e-01 6.971515e-02 2.675149e-02 1.060027e-01
## [2,] -2.302930e-01 4.240019e-01 3.235525e-02 7.621335e-02 1.937088e-01
## [3,] 6.971515e-02 3.235525e-02 3.338200e-01 -2.218598e-01 1.020704e-01
## [4,] 2.675149e-02 7.621335e-02 -2.218598e-01 4.079596e-01 1.029648e-01
## [5,] 1.060027e-01 1.937088e-01 1.020704e-01 1.029648e-01 2.997116e-01
## [6,] 9.646664e-02 1.085686e-01 1.119602e-01 1.860998e-01 2.050352e-01
## [7,] -3.221436e-02 -1.487315e-01 -9.972790e-02 -7.157314e-02 -1.809458e-01
## [8,] -1.212834e-01 -3.221494e-02 -8.531678e-02 -7.331967e-02 -1.534984e-01
## [9,] -8.459746e-02 -5.714882e-02 -5.305056e-02 -9.320243e-02 -1.417463e-01
## [10,] -1.549226e-01 -8.791262e-02 -1.683112e-01 -7.261134e-02 -2.428352e-01
## [11,] -1.460502e-01 -7.945741e-02 -7.871157e-02 -1.564158e-01 -2.255076e-01
## [12,] -7.019156e-02 -1.877674e-01 -1.450206e-01 -1.022485e-01 -2.579590e-01
## [13,] -1.612478e-01 -9.490402e-02 -1.542035e-01 -1.029546e-01 -2.561518e-01
## [14,] -6.971241e-02 -1.994368e-01 -9.510386e-02 -1.704180e-01 -2.691492e-01
## [15,] 1.658062e-17 -2.812342e-18 1.023140e-18 2.977664e-17 -2.177293e-17
## [16,] 6.812843e-18 4.417607e-17 -4.178285e-17 4.640101e-17 9.907072e-17
## [17,] 6.271884e-17 2.579012e-17 5.588952e-17 2.573204e-17 6.072600e-18
## [18,] -6.891580e-02 1.649668e-02 -6.483306e-02 1.183069e-02 -5.241912e-02
## [19,] -6.309608e-02 2.473172e-03 2.204704e-02 -9.483325e-02 -6.062291e-02
## [20,] 2.272884e-02 -8.840761e-02 -6.696106e-02 1.409512e-02 -6.567877e-02
## [21,] -6.599648e-02 1.972347e-02 -7.134430e-02 2.583526e-02 -4.627301e-02
## [22,] 3.242238e-02 -9.314286e-02 3.823423e-02 -9.978497e-02 -6.072048e-02
## [,6] [,7] [,8] [,9] [,10]
## [1,] 9.646664e-02 -3.221436e-02 -1.212834e-01 -8.459746e-02 -1.549226e-01
## [2,] 1.085686e-01 -1.487315e-01 -3.221494e-02 -5.714882e-02 -8.791262e-02
## [3,] 1.119602e-01 -9.972790e-02 -8.531678e-02 -5.305056e-02 -1.683112e-01
## [4,] 1.860998e-01 -7.157314e-02 -7.331967e-02 -9.320243e-02 -7.261134e-02
## [5,] 2.050352e-01 -1.809458e-01 -1.534984e-01 -1.417463e-01 -2.428352e-01
## [6,] 2.980600e-01 -1.713010e-01 -1.586365e-01 -1.462530e-01 -2.409226e-01
## [7,] -1.713010e-01 6.392168e-01 2.363105e-02 5.143787e-02 3.083020e-01
## [8,] -1.586365e-01 2.363105e-02 6.791768e-01 1.147783e-02 3.377197e-01
## [9,] -1.462530e-01 5.143787e-02 1.147783e-02 6.513700e-01 6.826402e-02
## [10,] -2.409226e-01 3.083020e-01 3.377197e-01 6.826402e-02 6.094300e-01
## [11,] -2.351274e-01 8.428748e-02 3.239027e-01 3.060956e-01 2.542925e-01
## [12,] -2.472691e-01 3.446969e-01 1.851742e-01 1.844147e-01 3.059343e-01
## [13,] -2.571581e-01 1.845548e-01 2.075013e-01 3.222297e-01 3.558104e-01
## [14,] -2.655219e-01 3.769941e-01 1.414741e-01 1.958175e-01 3.293019e-01
## [15,] -2.357044e-17 2.395790e-17 4.883271e-18 -2.708100e-17 1.803468e-17
## [16,] 1.286827e-16 -6.080268e-17 -1.220990e-16 -6.613368e-17 -1.055715e-16
## [17,] 1.048149e-17 -4.945417e-17 -3.400154e-17 -4.537775e-17 -6.495460e-17
## [18,] -5.300237e-02 -3.055726e-02 -3.616256e-02 6.671982e-02 -9.504654e-02
## [19,] -7.278621e-02 6.633679e-02 -3.255401e-02 -3.378277e-02 4.677847e-02
## [20,] -5.286593e-02 -2.046634e-02 1.368881e-02 6.777533e-03 2.956453e-02
## [21,] -4.550904e-02 2.655683e-02 2.393825e-02 -5.049508e-02 2.200864e-02
## [22,] -6.155074e-02 -4.187001e-02 3.108952e-02 1.078050e-02 -3.305105e-03
## [,11] [,12] [,13] [,14] [,15]
## [1,] -1.460502e-01 -7.019156e-02 -1.612478e-01 -6.971241e-02 -1.641941e-15
## [2,] -7.945741e-02 -1.877674e-01 -9.490402e-02 -1.994368e-01 -1.646013e-15
## [3,] -7.871157e-02 -1.450206e-01 -1.542035e-01 -9.510386e-02 4.117312e-16
## [4,] -1.564158e-01 -1.022485e-01 -1.029546e-01 -1.704180e-01 4.702551e-16
## [5,] -2.255076e-01 -2.579590e-01 -2.561518e-01 -2.691492e-01 1.566591e-15
## [6,] -2.351274e-01 -2.472691e-01 -2.571581e-01 -2.655219e-01 -5.109661e-16
## [7,] 8.428748e-02 3.446969e-01 1.845548e-01 3.769941e-01 1.598326e-16
## [8,] 3.239027e-01 1.851742e-01 2.075013e-01 1.414741e-01 3.106246e-17
## [9,] 3.060956e-01 1.844147e-01 3.222297e-01 1.958175e-01 -1.256668e-18
## [10,] 2.542925e-01 3.059343e-01 3.558104e-01 3.293019e-01 1.158805e-16
## [11,] 5.887043e-01 3.289698e-01 3.125844e-01 2.400521e-01 3.746834e-17
## [12,] 3.289698e-01 6.175962e-01 2.854138e-01 3.379905e-01 1.332151e-16
## [13,] 3.125844e-01 2.854138e-01 6.135308e-01 3.869710e-01 6.518561e-17
## [14,] 2.400521e-01 3.379905e-01 3.869710e-01 6.594414e-01 1.182485e-16
## [15,] -1.615002e-17 1.904230e-18 -1.674072e-18 1.360012e-17 4.285714e-01
## [16,] -1.424444e-16 -1.283217e-16 -9.498857e-17 -9.607831e-17 6.938894e-17
## [17,] -4.728673e-17 -5.206048e-17 -7.410980e-17 -6.773272e-17 2.949030e-17
## [18,] 4.080270e-02 2.505500e-02 2.703240e-02 9.434650e-03 -2.121626e-17
## [19,] -9.560950e-02 4.552758e-03 3.515618e-02 6.800707e-02 5.976514e-18
## [20,] -9.031318e-03 -9.165603e-02 4.421305e-02 3.536092e-02 -4.069808e-17
## [21,] 1.389874e-02 3.452602e-02 -9.191546e-02 -1.718182e-02 1.247080e-17
## [22,] 4.993938e-02 2.752225e-02 -1.448617e-02 -9.562082e-02 -3.932732e-17
## [,16] [,17] [,18] [,19] [,20]
## [1,] -2.325870e-16 -1.574980e-17 -6.891580e-02 -6.309608e-02 2.272884e-02
## [2,] -1.898534e-16 -4.427988e-17 1.649668e-02 2.473172e-03 -8.840761e-02
## [3,] 1.324522e-16 1.315660e-16 -6.483306e-02 2.204704e-02 -6.696106e-02
## [4,] 3.848965e-17 7.029055e-18 1.183069e-02 -9.483325e-02 1.409512e-02
## [5,] 2.572675e-16 3.622674e-17 -5.241912e-02 -6.062291e-02 -6.567877e-02
## [6,] -5.578429e-17 -7.509013e-17 -5.300237e-02 -7.278621e-02 -5.286593e-02
## [7,] -1.212212e-16 -1.196552e-17 -3.055726e-02 6.633679e-02 -2.046634e-02
## [8,] -5.859214e-17 -3.470884e-17 -3.616256e-02 -3.255401e-02 1.368881e-02
## [9,] 9.417605e-17 6.656983e-17 6.671982e-02 -3.378277e-02 6.777533e-03
## [10,] -9.870506e-17 -5.028129e-17 -9.504654e-02 4.677847e-02 2.956453e-02
## [11,] 3.671045e-18 3.963854e-17 4.080270e-02 -9.560950e-02 -9.031318e-03
## [12,] -8.981793e-17 2.475122e-17 2.505500e-02 4.552758e-03 -9.165603e-02
## [13,] 9.330102e-18 -5.458000e-18 2.703240e-02 3.515618e-02 4.421305e-02
## [14,] -5.351458e-17 -1.239952e-18 9.434650e-03 6.800707e-02 3.536092e-02
## [15,] -1.249001e-16 -2.602085e-17 -1.611237e-18 1.146659e-17 1.795250e-17
## [16,] 4.285714e-01 4.727121e-17 8.768364e-18 6.516845e-18 1.538432e-17
## [17,] -4.380177e-17 4.285714e-01 8.989137e-19 -1.359386e-18 3.590251e-18
## [18,] 3.104273e-17 6.352098e-18 3.298300e-01 1.866951e-02 2.391874e-02
## [19,] 2.529698e-18 -1.307603e-17 1.866951e-02 3.421286e-01 2.302693e-02
## [20,] 1.819295e-17 -4.029427e-17 2.391874e-02 2.302693e-02 3.362828e-01
## [21,] -4.051109e-17 -2.003227e-17 4.271666e-02 2.403364e-02 1.715788e-02
## [22,] 1.867247e-17 2.510164e-17 1.343650e-02 2.071270e-02 2.818508e-02
## [,21] [,22]
## [1,] -6.599648e-02 3.242238e-02
## [2,] 1.972347e-02 -9.314286e-02
## [3,] -7.134430e-02 3.823423e-02
## [4,] 2.583526e-02 -9.978497e-02
## [5,] -4.627301e-02 -6.072048e-02
## [6,] -4.550904e-02 -6.155074e-02
## [7,] 2.655683e-02 -4.187001e-02
## [8,] 2.393825e-02 3.108952e-02
## [9,] -5.049508e-02 1.078050e-02
## [10,] 2.200864e-02 -3.305105e-03
## [11,] 1.389874e-02 4.993938e-02
## [12,] 3.452602e-02 2.752225e-02
## [13,] -9.191546e-02 -1.448617e-02
## [14,] -1.718182e-02 -9.562082e-02
## [15,] 1.323016e-17 -1.324914e-18
## [16,] -9.098847e-19 -4.103621e-17
## [17,] 7.196728e-18 2.516784e-17
## [18,] 4.271666e-02 1.343650e-02
## [19,] 2.403364e-02 2.071270e-02
## [20,] 1.715788e-02 2.818508e-02
## [21,] 3.260662e-01 1.859699e-02
## [22,] 1.859699e-02 3.476402e-01
# solution
sol = gi_LHS %*% RHS
# print
sol
## [,1]
## [1,] 4.380219e+01
## [2,] 8.786764e+01
## [3,] 8.062658e+01
## [4,] 8.063976e+01
## [5,] 1.316698e+02
## [6,] 1.612663e+02
## [7,] 1.014757e+01
## [8,] -3.084153e+00
## [9,] -7.063418e+00
## [10,] 1.358074e+01
## [11,] -1.820697e+01
## [12,] -1.838678e+01
## [13,] 9.328432e+00
## [14,] 2.419362e+01
## [15,] -1.468082e-15
## [16,] 1.496534e-14
## [17,] -1.873529e-15
## [18,] 8.416979e+00
## [19,] -7.145576e+00
## [20,] -1.722850e+01
## [21,] -1.389647e+00
## [22,] 1.734674e+01
round(sol,3)
## [,1]
## [1,] 43.802
## [2,] 87.868
## [3,] 80.627
## [4,] 80.640
## [5,] 131.670
## [6,] 161.266
## [7,] 10.148
## [8,] -3.084
## [9,] -7.063
## [10,] 13.581
## [11,] -18.207
## [12,] -18.387
## [13,] 9.328
## [14,] 24.194
## [15,] 0.000
## [16,] 0.000
## [17,] 0.000
## [18,] 8.417
## [19,] -7.146
## [20,] -17.228
## [21,] -1.390
## [22,] 17.347
# Solution for fixed effects
no_of_levels_HYS = ncol(x1)
no_of_levels_parity = ncol(x2)
no_of_animals = ncol(Z)
# solution for HYS
sol_HYS =sol[1:no_of_levels_HYS]
sol_HYS
## [1] 43.80219 87.86764 80.62658 80.63976
# Solutions for the effect,parity
sol_parity =sol[(no_of_levels_HYS+1):(no_of_levels_HYS + no_of_levels_parity)]
sol_parity
## [1] 131.6698 161.2663
# solutions for the random animal effect
sol_animal = sol[(no_of_levels_HYS + no_of_levels_parity+1):(no_of_levels_HYS + no_of_levels_parity +
no_of_animals)]
sol_animal
## [1] 10.147570 -3.084153 -7.063418 13.580743 -18.206972 -18.386782 9.328432
## [8] 24.193618
# solutions for the permanent environment
sol_pe = sol[(no_of_levels_HYS + no_of_levels_parity + no_of_animals + 1):(
no_of_levels_HYS + no_of_levels_parity+no_of_animals*2)]
sol_pe
## [1] -1.468082e-15 1.496534e-14 -1.873529e-15 8.416979e+00 -7.145576e+00
## [6] -1.722850e+01 -1.389647e+00 1.734674e+01
# reliability(r2), accuracy(r), standard error of prediction(SEP)
# r2
# calculating the diagonal elements of the generalized inverse of the LHS
D = diag(gi_LHS[(( no_of_levels_HYS + no_of_levels_parity + 1):(
no_of_levels_HYS + no_of_levels_parity + no_of_animals)),((
no_of_levels_HYS + no_of_levels_parity + 1):(
no_of_levels_HYS + no_of_levels_parity + no_of_animals))])
D
## [1] 0.6392168 0.6791768 0.6513700 0.6094300 0.5887043 0.6175962 0.6135308
## [8] 0.6594414
r2 = 1 - D * alpha_1
r2
## [1] 0.10509649 0.04915244 0.08808198 0.14679806 0.17581391 0.13536537 0.14105693
## [8] 0.07678208
# r
r = sqrt(r2)
r
## [1] 0.3241859 0.2217035 0.2967861 0.3831423 0.4193017 0.3679203 0.3755755
## [8] 0.2770958
# SEP
sep = sqrt(D*sigma_2_e)
sep
## [1] 4.230611 4.360843 4.270639 4.130864 4.060015 4.158448 4.144739 4.297017
# Partitioning Breeding values
# yield deviation
yd1 = ginv( t(Z) %*% Z) %*% t(Z) %*% (y - x1 %*% sol_HYS - x2 %*% sol_parity - Z %*% sol_pe)
yd1
## [,1]
## [1,] 0.000000
## [2,] 0.000000
## [3,] 0.000000
## [4,] 23.400552
## [5,] -26.543478
## [6,] -38.486695
## [7,] 7.707178
## [8,] 44.431482
# numerator of n2
a2 = diag( t(Z) %*% Z)
a2
## [1] 0 0 0 2 2 2 2 2
# Parents average, progeny contribution
# parents avearge
pa1 = rep(0, no_of_animals)
# progeny contribution numerator
pc0 = rep(0, no_of_animals)
# numerator of n3, denominator of progeny contribution
a3 = rep(0, no_of_animals)
# numerator of n1
a1 = rep(0, no_of_animals)
# looping pedi
for (i in 1 : no_of_animals) {
sire = ped[i, 2]
dam = ped[i, 3]
if (sire == 0 & dam == 0) {
# both parents unknown
# PA
a1[i] = 1 * alpha_1
}
else if (sire != 0 & dam == 0) {
# PA
a1[i] = 4/3 * alpha_1
pa1[i] = sol_animal[sire]/2
# PC for sire
a3[sire] = a3[sire]+0.5 * alpha_1 * (2/3)
pc0[sire] = pc0[sire] + 0.5 * alpha_1 * (2/3) * (2 * sol_animal[i])
}
else if (sire == 0 & dam != 0) {
# dam known
# PA
a1[i] = 4/3 * alpha_1
pa1[i] = sol_animal[dam]/2
# PC for dam
a3[dam] = a3[dam] + 0.5 * alpha_1 * (2/3)
pc0[dam] = pc0[dam] + 0.5 * alpha_1 * (2/3) * (2 * sol_animal[i])
}
else {
# both parents known
# PA
a1[i] = 2 * alpha_1
pa1[i] = (sol_animal[sire] + sol_animal[dam])/2
# PC for sire
a3[sire] = a3[sire] + 0.5 * alpha_1
pc0[sire] = pc0[sire] + 0.5 * alpha_1 * (2 * sol_animal[i] - sol_animal[dam])
# PC for dam
a3[dam] = a3[dam] + 0.5 * alpha_1
pc0[dam] = pc0[dam] + 0.5 * alpha_1 * (2 * sol_animal[i] - sol_animal[sire])
}
}
a1
## [1] 1.4 1.4 1.4 2.8 2.8 2.8 2.8 2.8
pa1
## [1] 0.000000 0.000000 0.000000 3.531709 -5.073785 -4.029701 3.258663
## [8] 9.738001
a3
## [1] 2.1 1.4 1.4 0.7 0.7 0.0 0.7 0.0
pc0
## [1] 35.516497 -8.635628 -19.777569 18.004198 -32.844794 0.000000 26.767766
## [8] 0.000000
# denominator of n1, n2, n3, diagonal of animals in LHS
n_de = a1 + a2 + a3
n_de
## [1] 3.5 2.8 2.8 5.5 5.5 4.8 5.5 4.8
# parents average fraction of breeding values
pa = pa1 * (a1 / n_de)
pa
## [1] 0.000000 0.000000 0.000000 1.797961 -2.583018 -2.350659 1.658956
## [8] 5.680501
# yield deviation fraction of breeding values
yd = yd1 * (a2 / n_de)
yd
## [,1]
## [1,] 0.000000
## [2,] 0.000000
## [3,] 0.000000
## [4,] 8.509292
## [5,] -9.652174
## [6,] -16.036123
## [7,] 2.802610
## [8,] 18.513118
# progeny contribution
pc1 = pc0 / a3
pc1[is.nan(pc1) == TRUE] = 0
pc1
## [1] 16.912617 -6.168306 -14.126835 25.720282 -46.921134 0.000000 38.239666
## [8] 0.000000
# progeny contribution fraction of breeding values
pc = pc1 * (a3 / n_de)
pc
## [1] 10.147570 -3.084153 -7.063418 3.273490 -5.971781 0.000000 4.866867
## [8] 0.000000
# Progeny(Daughter) Yield Deviation(PYD, DYD)
# n2 of progeny
n2prog = a2 / (a1 + a2)
n2prog
## [1] 0.0000000 0.0000000 0.0000000 0.4166667 0.4166667 0.4166667 0.4166667
## [8] 0.4166667
# numerator of dyd : summation of u of progeny * n2 of progeny * (2 * YD - bv of mate)
dyd_n = rep(0, no_of_animals)
dyd_n
## [1] 0 0 0 0 0 0 0 0
# denominator of dyd : summation of u of progeny * n2 of progeny
dyd_d = rep(0, no_of_animals)
dyd_d
## [1] 0 0 0 0 0 0 0 0
# looping pedi
for (i in 1 : no_of_animals) {
sire = ped[i, 2]
dam = ped[i, 3]
if (sire == 0 & dam == 0) {
# both parents unknown
} else if (sire != 0 & dam == 0) {
# dyd_n
dyd_n[sire] = dyd_n[sire] + n2prog[i] * 2 / 3 * 2 * yd1[i]
# dyd_d
dyd_d[sire] = dyd_d[sire] + n2prog[i] * 2 / 3
} else if (sire == 0 & dam != 0) {
# dam known
# dyd_n
dyd_n[dam] = dyd_n[dam] + n2prog[i] * 2 / 3 * 2 * yd1[i]
# dyd_d
dyd_d[dam] = dyd_d[dam]
n2prog[i] * 2 / 3
} else {
# both parents known
# dyd_n
dyd_n[sire] = dyd_n[sire] + n2prog[i] * (2 * yd1[i] - sol_animal[dam])
dyd_n[dam] = dyd_n[dam] + n2prog[i] * (2 * yd1[i] - sol_animal[sire])
# dyd_d
dyd_d[sire] = dyd_d[sire] + n2prog[i]
dyd_d[dam] = dyd_d[dam] + n2prog[i]
}
}
dyd_n
## [1] 29.438905 -3.904169 -20.070496 9.365739 -36.300400 0.000000 32.798081
## [8] 0.000000
dyd_d
## [1] 1.2500000 0.8333333 0.8333333 0.4166667 0.4166667 0.0000000 0.4166667
## [8] 0.0000000
# dyd
dyd = dyd_n / dyd_d
dyd[is.nan(dyd) == TRUE] = 0
dyd
## [1] 23.551124 -4.685002 -24.084595 22.477774 -87.120960 0.000000 78.715394
## [8] 0.000000
# breeding values and fractions
Repeatability_Model_results <- data.frame(animal = ped[,1], animal_bv = sol_animal,
rel = r2, acc = r, sep = sep, pa = pa, yd = yd,
pc = pc, sum_of_fr = pa + yd + pc, dyd = dyd)
Repeatability_Model_results
## animal animal_bv rel acc sep pa yd
## 1 1 10.147570 0.10509649 0.3241859 4.230611 0.000000 0.000000
## 2 2 -3.084153 0.04915244 0.2217035 4.360843 0.000000 0.000000
## 3 3 -7.063418 0.08808198 0.2967861 4.270639 0.000000 0.000000
## 4 4 13.580743 0.14679806 0.3831423 4.130864 1.797961 8.509292
## 5 5 -18.206972 0.17581391 0.4193017 4.060015 -2.583018 -9.652174
## 6 6 -18.386782 0.13536537 0.3679203 4.158448 -2.350659 -16.036123
## 7 7 9.328432 0.14105693 0.3755755 4.144739 1.658956 2.802610
## 8 8 24.193618 0.07678208 0.2770958 4.297017 5.680501 18.513118
## pc sum_of_fr dyd
## 1 10.147570 10.147570 23.551124
## 2 -3.084153 -3.084153 -4.685002
## 3 -7.063418 -7.063418 -24.084595
## 4 3.273490 13.580743 22.477774
## 5 -5.971781 -18.206972 -87.120960
## 6 0.000000 -18.386782 0.000000
## 7 4.866867 9.328432 78.715394
## 8 0.000000 24.193618 0.000000