Multiple trait animal model, same model and no missing record
Linear Models for the Prediction of Animal Breeding Values, 3rd Edtion
Raphael Mrode
Example 5.1 p72
자료
4 1 4.5 6.8
5 2 2.9 5.0
6 2 3.9 6.8
7 1 3.5 6.0
8 1 5.0 7.5animal, sex, wwg(pre-weaning gain), pwg(post-weaning gain)
mt01_data.txt로 저장
혈통
1 0 0
2 0 0
3 0 0
4 1 0
5 3 2
6 1 2
7 4 5
8 3 6
자료와 혈통을 읽어 육종가를 구하는 R 프로그램이다.
신뢰도, 정확도, Standard error of prediction을 구한다.
육종가를 PA(parenta average), YD(yield devation), PC(progeny contribution)으로 구한다.
개체의 DYD(daughter yield deviation)을 구한다.
# multiple trait animal model, same model and no missing record - Date : 2020-07-03
# Linear Models for the Prediction of Animal Breeding Values, 3rd Edtion
# Raphael Mrode
# Example 5.1 p72
# MASS packages
if (!("MASS" %in% installed.packages())){
install.packages('MASS', dependencies = TRUE)
}
library(MASS)
# Matrix packages
if (!("Matrix" %in% installed.packages())){
install.packages('Matrix', dependencies = TRUE)
}
library(Matrix)
# functions
# 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)
# 각 수준의 인덱스 값을 저장
va = as.numeric(vf)
mrow = length(va)
mcol = length(levels(vf))
}
# 0으로 채워진 X 준비
X = matrix(data = c(0), nrow = mrow, ncol = mcol)
for (i in 1:mrow) {
ic = va[i]
X[i, ic] = 1
}
return(X)
}
# function to make inverse of numerator relationship matrix
ainv = function(pedi) {
n = nrow(pedi)
Ainv = matrix(c(0), nrow = n, ncol = n)
for (i in 1:n) {
animal = pedi[i, 1]
sire = pedi[i, 2]
dam = pedi[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)
}
# set working directory
setwd("D:/temp/04_multiple_traits_01_R")
# print working directory
getwd()
# no_of_trts
no_of_trts = 2
# prdigree and data
# read pedigree : animal sire dam
pedi = read.table("mt01_pedi.txt", header = FALSE, sep = "", col.names = c("animal", "sire", "dam"),stringsAsFactors = FALSE)
pedi = pedi[order(pedi$animal),]
# print
pedi
# read data : animal, dam, sex, weaning_weight
data = read.table("mt01_data.txt", header = FALSE, sep = "", col.names = c("animal", "sex", "wwg", "pwg"), stringsAsFactors = FALSE)
# print
data
# number of data
no_data = nrow(data)
# print
no_data
# variance component and ratio
G = matrix(c(20, 18, 18, 40), 2, 2)
R = matrix(c(40, 11, 11, 30), 2, 2)
# print
G
R
# inverse of G and R
Gi = ginv(G)
Ri = ginv(R)
# print
Gi
Ri
# full matrix of R inverse
Rif = kronecker(Ri, diag(no_data))
# print
Rif
# design matrix
# design matrix of 1st fixed effect
X1 = desgn(data[, 2])
# block diagonal matirx
X = as.matrix(bdiag(X1, X1))
# print
X
# number of levels of 1st fixed effect
no_lvl_f1 = ncol(X1)
# print
no_lvl_f1
# desing matrix of animal effect
Z1 = desgn(data[, 1])
# block diagonal matirx
Z = as.matrix(bdiag(Z1, Z1))
# print
Z
# number of animals
no_animals = ncol(Z1)
# print
no_animals
# observation
y1 = as.matrix(data[, 3])
y2 = as.matrix(data[, 4])
y = rbind(y1, y2)
# print
y1
y2
y
# inverse matrix of NRM
ai = ainv(pedi)
# print
ai
# LHS, RHS
# LHS construction
XPX = t(X) %*% Rif %*% X
XPZ = t(X) %*% Rif %*% Z
ZPX = t(Z) %*% Rif %*% X
ZPZ = t(Z) %*% Rif %*% Z
#print
XPX
XPZ
ZPX
ZPZ
LHS = rbind(
cbind(XPX, XPZ),
cbind(ZPX, ZPZ + kronecker(Gi, ai))
)
# print
LHS
# RHS construction
XPy = t(X) %*% Rif %*% y
ZPy = t(Z) %*% Rif %*% y
# print
XPy
ZPy
RHS = rbind(XPy,
ZPy
)
# print
RHS
# Solutions
# generalized inverse of LHS
gi_LHS = ginv(LHS)
# print
gi_LHS
# solution
sol = gi_LHS %*% RHS
# print
sol
# solutions for fixed effects and animal effects
sol_t1_f1 = as.matrix(sol[1 : no_lvl_f1])
sol_t2_f1 = as.matrix(sol[(no_lvl_f1 + 1) : (no_lvl_f1 * 2)])
sol_f1 = cbind(sol_t1_f1, sol_t2_f1)
#print
sol_f1
# breedinv value
sol_t1_bv = sol[(no_lvl_f1 * 2 + 1) : (no_lvl_f1 * 2 + no_animals)]
sol_t2_bv = sol[(no_lvl_f1 * 2 + no_animals + 1) : (no_lvl_f1 * 2 + no_animals *2)]
sol_bv = cbind(sol_t1_bv, sol_t2_bv)
# print
sol_bv
# reliability(r2), accuracy(r), standard error of prediction(SEP)
# diagonal elements of the generalized inverse of LHS for animal equation
d_ani_t1 = diag(gi_LHS[(no_lvl_f1 * 2 + 1) : (no_lvl_f1 * 2 + no_animals),
(no_lvl_f1 * 2 + 1) : (no_lvl_f1 * 2 + no_animals)])
d_ani_t2 = diag(gi_LHS[(no_lvl_f1 * 2 + no_animals + 1) : (no_lvl_f1 * 2 + no_animals * 2),
(no_lvl_f1 * 2 + no_animals + 1) : (no_lvl_f1 * 2 + no_animals * 2)])
d_ani = cbind(d_ani_t1, d_ani_t2)
# print
d_ani
# reliability
rel = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
for(i in 1 : no_animals){
rel[i, ] = 1 - d_ani[i,] / diag(G)
}
# print
rel
# accuracy
acc = sqrt(rel)
# print
acc
# standard error of prediction(SEP)
sep = sqrt(d_ani)
# 확인
sep
# partitioning of breeding values
# yield deviation
yd1 = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
# numerator of n2
a2 = array(rep(0, no_of_trts * no_of_trts * no_animals), dim = c(no_of_trts, no_of_trts, no_animals))
# looping data
for (i in 1 : no_data) {
yd1[data[i, 1], ] = as.matrix(data[i, c(3,4)] - sol_f1[data[i,2],])
a2[,,data[i, 1]] = Ri
}
# print
yd1
a2
# Parents average, progeny contribution
# parents avearge
pa1 = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
# progeny contribution numerator
pc0 = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
# numerator of n3, denominator of progeny contribution
a3 = array(rep(0, no_of_trts * no_of_trts * no_animals), dim = c(no_of_trts, no_of_trts, no_animals))
# numerator of n1
a1 = array(rep(0, no_of_trts * no_of_trts * no_animals), dim = c(no_of_trts, no_of_trts, no_animals))
# looping pedi
for (i in 1 : no_animals) {
sire = pedi[i, 2]
dam = pedi[i, 3]
if (sire == 0 & dam == 0) {
# both parents unknown
# PA
a1[ , , i] = 1 * Gi
} else if (sire != 0 & dam == 0) {
# 개체의 부만 알고 있을 경우
# PA
a1[ , , i] = 4/3 * Gi
pa1[i, ] = sol_bv[sire, ]/2
# PC for sire
a3[ , , sire] = a3[ , , sire] + 0.5 * Gi * (2/3)
pc0[sire, ] = pc0[sire, ] + (0.5 * Gi * (2/3)) %*% (2 * sol_bv[i,])
} else if (sire == 0 & dam != 0) {
# dam known
# PA
a1[ , , i] = 4/3 * Gi
pa1[i, ] = sol_bv[dam, ]/2
# PC for dam
a3[ , , dam] = a3[ , , dam] + 0.5 * Gi * (2/3)
pc0[dam, ] = pc0[dam, ] + (0.5 * Gi * (2/3)) %*% (2 * sol_bv[i])
} else {
# both parents known
# PA
a1[ , , i] = 2 * Gi
pa1[i, ] = (sol_bv[sire, ] + sol_bv[dam,])/2
# PC for sire
a3[ , , sire] = a3[ , , sire] + 0.5 * Gi
pc0[sire, ] = pc0[sire, ] + (0.5 * Gi) %*% (2 * sol_bv[i, ] - sol_bv[dam, ])
# PC for dam
a3[ , , dam] = a3[ , , dam] + 0.5 * Gi
pc0[dam, ] = pc0[dam, ] + (0.5 * Gi) %*% (2 * sol_bv[i, ] - sol_bv[sire, ])
}
}
# print
a1
pa1
a3
pc0
# denominator of n1, n2, n3, diagonal of animals in LHS
n_de = a1 + a2 + a3
# print
n_de
# parents average fraction of breeding values
pa = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
for (i in 1 : no_animals) {
pa[i, ] = ginv(n_de[ , , i]) %*% a1[, , i] %*% pa1[i, ]
}
# print
pa
# yield deviation fraction of breeding values
yd = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
for (i in 1 : no_animals) {
yd[i, ] = ginv(n_de[ , , i]) %*% a2[, , i] %*% yd1[i, ]
}
# print
yd
# progeny contribution
pc1 = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
for (i in 1 : no_animals) {
pc1[i, ] = ginv(a3[ , , i]) %*% pc0[i, ]
}
pc1[is.nan(pc1) == TRUE] = 0
pc1
# progeny contribution fraction of breeding values
pc = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
for (i in 1 : no_animals) {
pc[i, ] = ginv(n_de[ , , i]) %*% a3[, , i] %*% pc1[i, ]
}
# print
pc
# Progeny(Daughter) Yield Deviation(PYD, DYD)
# n2 of progeny
n2prog = array(rep(0, no_of_trts * no_of_trts * no_animals), dim = c(no_of_trts, no_of_trts, no_animals))
for (i in 1 : no_animals) {
n2prog[ , , i] = ginv((a1 + a2)[, , i]) %*% a2[, , i]
}
# print
n2prog
# numerator of dyd : summation of u of progeny * n2 of progeny * (2 * YD - bv of mate)
dyd_n = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
# denominator of dyd : summation of u of progeny * n2 of progeny
dyd_d = array(rep(0, no_of_trts * no_of_trts * no_animals), dim = c(no_of_trts, no_of_trts, no_animals))
# looping pedi
for (i in 1 : no_animals) {
# i = 5
sire = pedi[i, 2]
dam = pedi[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_bv[dam, ])
dyd_n[dam, ] = dyd_n[dam, ] + n2prog[ , , i] %*% (2 * yd1[i, ] - sol_bv[sire, ])
# dyd_d
dyd_d[ , , sire] = dyd_d[ , , sire] + n2prog[ , , i]
dyd_d[ , , dam] = dyd_d[ , , dam] + n2prog[ , , i]
}
}
# print
dyd_n
dyd_d
# dyd
dyd = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
# looping pedi
for (i in 1 : no_animals) {
dyd[i, ] = ginv(dyd_d[ , , i]) %*% dyd_n[i, ]
}
dyd[is.nan(dyd) == TRUE] = 0
# print
dyd
# breeding values and fractions
mt1_result = data.frame(animal = pedi[,1], animal_bv = sol_bv, rel = rel, acc = acc, sep = sep, pa = pa, yd = yd, pc = pc, sum_of_fr = pa + yd + pc, dyd = dyd)
# print
mt1_result
# 파일로 출력, 분리자는 ",", 따옴표 없고
output_filename = gsub("[ -]", "", paste("mt1_result_", Sys.Date(), ".csv"))
write.table(mt1_result, output_filename, sep=",", row.names = FALSE, quote = FALSE)
다음은 R 프로그램을 실행한 결과이다.
> # multiple trait animal model, same model and no missing record - Date : 2020-07-03
>
> # Linear Models for the Prediction of Animal Breeding Values, 3rd Edtion
>
> # Raphael Mrode
>
> # Example 5.1 p72
>
> # MASS packages
> if (!("MASS" %in% installed.packages())){
+ install.packages('MASS', dependencies = TRUE)
+ }
> library(MASS)
>
> # Matrix packages
> if (!("Matrix" %in% installed.packages())){
+ install.packages('Matrix', dependencies = TRUE)
+ }
> library(Matrix)
>
> # functions
>
> # 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)
+ # 각 수준의 인덱스 값을 저장
+ va = as.numeric(vf)
+ mrow = length(va)
+ mcol = length(levels(vf))
+ }
+
+ # 0으로 채워진 X 준비
+ X = matrix(data = c(0), nrow = mrow, ncol = mcol)
+
+ for (i in 1:mrow) {
+ ic = va[i]
+ X[i, ic] = 1
+ }
+ return(X)
+ }
>
> # function to make inverse of numerator relationship matrix
> ainv = function(pedi) {
+ n = nrow(pedi)
+ Ainv = matrix(c(0), nrow = n, ncol = n)
+
+ for (i in 1:n) {
+ animal = pedi[i, 1]
+ sire = pedi[i, 2]
+ dam = pedi[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)
+ }
>
> # set working directory
> setwd("D:/temp/04_multiple_traits_01_R")
>
> # print working directory
> getwd()
[1] "D:/temp/04_multiple_traits_01_R"
>
> # no_of_trts
> no_of_trts = 2
>
> # prdigree and data
>
> # read pedigree : animal sire dam
> pedi = read.table("mt01_pedi.txt", header = FALSE, sep = "", col.names = c("animal", "sire", "dam"),stringsAsFactors = FALSE)
> pedi = pedi[order(pedi$animal),]
>
> # print
> pedi
animal sire dam
1 1 0 0
2 2 0 0
3 3 0 0
4 4 1 0
5 5 3 2
6 6 1 2
7 7 4 5
8 8 3 6
>
> # read data : animal, dam, sex, weaning_weight
> data = read.table("mt01_data.txt", header = FALSE, sep = "", col.names = c("animal", "sex", "wwg", "pwg"), stringsAsFactors = FALSE)
>
> # print
> data
animal sex wwg pwg
1 4 1 4.5 6.8
2 5 2 2.9 5.0
3 6 2 3.9 6.8
4 7 1 3.5 6.0
5 8 1 5.0 7.5
>
> # number of data
> no_data = nrow(data)
>
> # print
> no_data
[1] 5
>
> # variance component and ratio
> G = matrix(c(20, 18, 18, 40), 2, 2)
> R = matrix(c(40, 11, 11, 30), 2, 2)
>
> # print
> G
[,1] [,2]
[1,] 20 18
[2,] 18 40
> R
[,1] [,2]
[1,] 40 11
[2,] 11 30
>
> # inverse of G and R
> Gi = ginv(G)
> Ri = ginv(R)
>
> # print
> Gi
[,1] [,2]
[1,] 0.08403361 -0.03781513
[2,] -0.03781513 0.04201681
> Ri
[,1] [,2]
[1,] 0.02780352 -0.01019462
[2,] -0.01019462 0.03707136
>
> # full matrix of R inverse
> Rif = kronecker(Ri, diag(no_data))
>
> # print
> Rif
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 0.02780352 0.00000000 0.00000000 0.00000000 0.00000000 -0.01019462 0.00000000
[2,] 0.00000000 0.02780352 0.00000000 0.00000000 0.00000000 0.00000000 -0.01019462
[3,] 0.00000000 0.00000000 0.02780352 0.00000000 0.00000000 0.00000000 0.00000000
[4,] 0.00000000 0.00000000 0.00000000 0.02780352 0.00000000 0.00000000 0.00000000
[5,] 0.00000000 0.00000000 0.00000000 0.00000000 0.02780352 0.00000000 0.00000000
[6,] -0.01019462 0.00000000 0.00000000 0.00000000 0.00000000 0.03707136 0.00000000
[7,] 0.00000000 -0.01019462 0.00000000 0.00000000 0.00000000 0.00000000 0.03707136
[8,] 0.00000000 0.00000000 -0.01019462 0.00000000 0.00000000 0.00000000 0.00000000
[9,] 0.00000000 0.00000000 0.00000000 -0.01019462 0.00000000 0.00000000 0.00000000
[10,] 0.00000000 0.00000000 0.00000000 0.00000000 -0.01019462 0.00000000 0.00000000
[,8] [,9] [,10]
[1,] 0.00000000 0.00000000 0.00000000
[2,] 0.00000000 0.00000000 0.00000000
[3,] -0.01019462 0.00000000 0.00000000
[4,] 0.00000000 -0.01019462 0.00000000
[5,] 0.00000000 0.00000000 -0.01019462
[6,] 0.00000000 0.00000000 0.00000000
[7,] 0.00000000 0.00000000 0.00000000
[8,] 0.03707136 0.00000000 0.00000000
[9,] 0.00000000 0.03707136 0.00000000
[10,] 0.00000000 0.00000000 0.03707136
>
> # design matrix
>
> # design matrix of 1st fixed effect
> X1 = desgn(data[, 2])
>
> # block diagonal matirx
> X = as.matrix(bdiag(X1, X1))
>
> # print
> X
[,1] [,2] [,3] [,4]
[1,] 1 0 0 0
[2,] 0 1 0 0
[3,] 0 1 0 0
[4,] 1 0 0 0
[5,] 1 0 0 0
[6,] 0 0 1 0
[7,] 0 0 0 1
[8,] 0 0 0 1
[9,] 0 0 1 0
[10,] 0 0 1 0
>
> # number of levels of 1st fixed effect
> no_lvl_f1 = ncol(X1)
>
> # print
> no_lvl_f1
[1] 2
>
> # desing matrix of animal effect
> Z1 = desgn(data[, 1])
>
> # block diagonal matirx
> Z = as.matrix(bdiag(Z1, Z1))
>
> # print
> Z
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
[1,] 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
[2,] 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
[3,] 0 0 0 0 0 1 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
[5,] 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
[6,] 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
[7,] 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
[8,] 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
[9,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
[10,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
>
> # number of animals
> no_animals = ncol(Z1)
>
> # print
> no_animals
[1] 8
>
> # observation
> y1 = as.matrix(data[, 3])
> y2 = as.matrix(data[, 4])
> y = rbind(y1, y2)
>
> # print
> y1
[,1]
[1,] 4.5
[2,] 2.9
[3,] 3.9
[4,] 3.5
[5,] 5.0
> y2
[,1]
[1,] 6.8
[2,] 5.0
[3,] 6.8
[4,] 6.0
[5,] 7.5
> y
[,1]
[1,] 4.5
[2,] 2.9
[3,] 3.9
[4,] 3.5
[5,] 5.0
[6,] 6.8
[7,] 5.0
[8,] 6.8
[9,] 6.0
[10,] 7.5
>
> # inverse matrix of NRM
> ai = ainv(pedi)
>
> # print
> ai
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1.8333333 0.5 0.0 -0.6666667 0.0 -1.0 0 0
[2,] 0.5000000 2.0 0.5 0.0000000 -1.0 -1.0 0 0
[3,] 0.0000000 0.5 2.0 0.0000000 -1.0 0.5 0 -1
[4,] -0.6666667 0.0 0.0 1.8333333 0.5 0.0 -1 0
[5,] 0.0000000 -1.0 -1.0 0.5000000 2.5 0.0 -1 0
[6,] -1.0000000 -1.0 0.5 0.0000000 0.0 2.5 0 -1
[7,] 0.0000000 0.0 0.0 -1.0000000 -1.0 0.0 2 0
[8,] 0.0000000 0.0 -1.0 0.0000000 0.0 -1.0 0 2
>
> # LHS, RHS
>
> # LHS construction
> XPX = t(X) %*% Rif %*% X
> XPZ = t(X) %*% Rif %*% Z
> ZPX = t(Z) %*% Rif %*% X
> ZPZ = t(Z) %*% Rif %*% Z
>
> #print
> XPX
[,1] [,2] [,3] [,4]
[1,] 0.08341057 0.00000000 -0.03058387 0.00000000
[2,] 0.00000000 0.05560704 0.00000000 -0.02038925
[3,] -0.03058387 0.00000000 0.11121409 0.00000000
[4,] 0.00000000 -0.02038925 0.00000000 0.07414272
> XPZ
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0 0 0 0.02780352 0.00000000 0.00000000 0.02780352 0.02780352 0 0
[2,] 0 0 0 0.00000000 0.02780352 0.02780352 0.00000000 0.00000000 0 0
[3,] 0 0 0 -0.01019462 0.00000000 0.00000000 -0.01019462 -0.01019462 0 0
[4,] 0 0 0 0.00000000 -0.01019462 -0.01019462 0.00000000 0.00000000 0 0
[,11] [,12] [,13] [,14] [,15] [,16]
[1,] 0 -0.01019462 0.00000000 0.00000000 -0.01019462 -0.01019462
[2,] 0 0.00000000 -0.01019462 -0.01019462 0.00000000 0.00000000
[3,] 0 0.03707136 0.00000000 0.00000000 0.03707136 0.03707136
[4,] 0 0.00000000 0.03707136 0.03707136 0.00000000 0.00000000
> ZPX
[,1] [,2] [,3] [,4]
[1,] 0.00000000 0.00000000 0.00000000 0.00000000
[2,] 0.00000000 0.00000000 0.00000000 0.00000000
[3,] 0.00000000 0.00000000 0.00000000 0.00000000
[4,] 0.02780352 0.00000000 -0.01019462 0.00000000
[5,] 0.00000000 0.02780352 0.00000000 -0.01019462
[6,] 0.00000000 0.02780352 0.00000000 -0.01019462
[7,] 0.02780352 0.00000000 -0.01019462 0.00000000
[8,] 0.02780352 0.00000000 -0.01019462 0.00000000
[9,] 0.00000000 0.00000000 0.00000000 0.00000000
[10,] 0.00000000 0.00000000 0.00000000 0.00000000
[11,] 0.00000000 0.00000000 0.00000000 0.00000000
[12,] -0.01019462 0.00000000 0.03707136 0.00000000
[13,] 0.00000000 -0.01019462 0.00000000 0.03707136
[14,] 0.00000000 -0.01019462 0.00000000 0.03707136
[15,] -0.01019462 0.00000000 0.03707136 0.00000000
[16,] -0.01019462 0.00000000 0.03707136 0.00000000
> ZPZ
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0 0 0 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0 0
[2,] 0 0 0 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0 0
[3,] 0 0 0 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0 0
[4,] 0 0 0 0.02780352 0.00000000 0.00000000 0.00000000 0.00000000 0 0
[5,] 0 0 0 0.00000000 0.02780352 0.00000000 0.00000000 0.00000000 0 0
[6,] 0 0 0 0.00000000 0.00000000 0.02780352 0.00000000 0.00000000 0 0
[7,] 0 0 0 0.00000000 0.00000000 0.00000000 0.02780352 0.00000000 0 0
[8,] 0 0 0 0.00000000 0.00000000 0.00000000 0.00000000 0.02780352 0 0
[9,] 0 0 0 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0 0
[10,] 0 0 0 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0 0
[11,] 0 0 0 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0 0
[12,] 0 0 0 -0.01019462 0.00000000 0.00000000 0.00000000 0.00000000 0 0
[13,] 0 0 0 0.00000000 -0.01019462 0.00000000 0.00000000 0.00000000 0 0
[14,] 0 0 0 0.00000000 0.00000000 -0.01019462 0.00000000 0.00000000 0 0
[15,] 0 0 0 0.00000000 0.00000000 0.00000000 -0.01019462 0.00000000 0 0
[16,] 0 0 0 0.00000000 0.00000000 0.00000000 0.00000000 -0.01019462 0 0
[,11] [,12] [,13] [,14] [,15] [,16]
[1,] 0 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
[2,] 0 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
[3,] 0 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
[4,] 0 -0.01019462 0.00000000 0.00000000 0.00000000 0.00000000
[5,] 0 0.00000000 -0.01019462 0.00000000 0.00000000 0.00000000
[6,] 0 0.00000000 0.00000000 -0.01019462 0.00000000 0.00000000
[7,] 0 0.00000000 0.00000000 0.00000000 -0.01019462 0.00000000
[8,] 0 0.00000000 0.00000000 0.00000000 0.00000000 -0.01019462
[9,] 0 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
[10,] 0 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
[11,] 0 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
[12,] 0 0.03707136 0.00000000 0.00000000 0.00000000 0.00000000
[13,] 0 0.00000000 0.03707136 0.00000000 0.00000000 0.00000000
[14,] 0 0.00000000 0.00000000 0.03707136 0.00000000 0.00000000
[15,] 0 0.00000000 0.00000000 0.00000000 0.03707136 0.00000000
[16,] 0 0.00000000 0.00000000 0.00000000 0.00000000 0.03707136
>
> LHS = rbind(
+ cbind(XPX, XPZ),
+ cbind(ZPX, ZPZ + kronecker(Gi, ai))
+ )
>
> # print
> LHS
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 0.08341057 0.00000000 -0.03058387 0.00000000 0.00000000 0.00000000 0.00000000
[2,] 0.00000000 0.05560704 0.00000000 -0.02038925 0.00000000 0.00000000 0.00000000
[3,] -0.03058387 0.00000000 0.11121409 0.00000000 0.00000000 0.00000000 0.00000000
[4,] 0.00000000 -0.02038925 0.00000000 0.07414272 0.00000000 0.00000000 0.00000000
[5,] 0.00000000 0.00000000 0.00000000 0.00000000 0.15406162 0.04201681 0.00000000
[6,] 0.00000000 0.00000000 0.00000000 0.00000000 0.04201681 0.16806723 0.04201681
[7,] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.04201681 0.16806723
[8,] 0.02780352 0.00000000 -0.01019462 0.00000000 -0.05602241 0.00000000 0.00000000
[9,] 0.00000000 0.02780352 0.00000000 -0.01019462 0.00000000 -0.08403361 -0.08403361
[10,] 0.00000000 0.02780352 0.00000000 -0.01019462 -0.08403361 -0.08403361 0.04201681
[11,] 0.02780352 0.00000000 -0.01019462 0.00000000 0.00000000 0.00000000 0.00000000
[12,] 0.02780352 0.00000000 -0.01019462 0.00000000 0.00000000 0.00000000 -0.08403361
[13,] 0.00000000 0.00000000 0.00000000 0.00000000 -0.06932773 -0.01890756 0.00000000
[14,] 0.00000000 0.00000000 0.00000000 0.00000000 -0.01890756 -0.07563025 -0.01890756
[15,] 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 -0.01890756 -0.07563025
[16,] -0.01019462 0.00000000 0.03707136 0.00000000 0.02521008 0.00000000 0.00000000
[17,] 0.00000000 -0.01019462 0.00000000 0.03707136 0.00000000 0.03781513 0.03781513
[18,] 0.00000000 -0.01019462 0.00000000 0.03707136 0.03781513 0.03781513 -0.01890756
[19,] -0.01019462 0.00000000 0.03707136 0.00000000 0.00000000 0.00000000 0.00000000
[20,] -0.01019462 0.00000000 0.03707136 0.00000000 0.00000000 0.00000000 0.03781513
[,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 0.02780352 0.00000000 0.00000000 0.02780352 0.02780352 0.00000000 0.00000000
[2,] 0.00000000 0.02780352 0.02780352 0.00000000 0.00000000 0.00000000 0.00000000
[3,] -0.01019462 0.00000000 0.00000000 -0.01019462 -0.01019462 0.00000000 0.00000000
[4,] 0.00000000 -0.01019462 -0.01019462 0.00000000 0.00000000 0.00000000 0.00000000
[5,] -0.05602241 0.00000000 -0.08403361 0.00000000 0.00000000 -0.06932773 -0.01890756
[6,] 0.00000000 -0.08403361 -0.08403361 0.00000000 0.00000000 -0.01890756 -0.07563025
[7,] 0.00000000 -0.08403361 0.04201681 0.00000000 -0.08403361 0.00000000 -0.01890756
[8,] 0.18186515 0.04201681 0.00000000 -0.08403361 0.00000000 0.02521008 0.00000000
[9,] 0.04201681 0.23788756 0.00000000 -0.08403361 0.00000000 0.00000000 0.03781513
[10,] 0.00000000 0.00000000 0.23788756 0.00000000 -0.08403361 0.03781513 0.03781513
[11,] -0.08403361 -0.08403361 0.00000000 0.19587075 0.00000000 0.00000000 0.00000000
[12,] 0.00000000 0.00000000 -0.08403361 0.00000000 0.19587075 0.00000000 0.00000000
[13,] 0.02521008 0.00000000 0.03781513 0.00000000 0.00000000 0.07703081 0.02100840
[14,] 0.00000000 0.03781513 0.03781513 0.00000000 0.00000000 0.02100840 0.08403361
[15,] 0.00000000 0.03781513 -0.01890756 0.00000000 0.03781513 0.00000000 0.02100840
[16,] -0.07952236 -0.01890756 0.00000000 0.03781513 0.00000000 -0.02801120 0.00000000
[17,] -0.01890756 -0.10473244 0.00000000 0.03781513 0.00000000 0.00000000 -0.04201681
[18,] 0.00000000 0.00000000 -0.10473244 0.00000000 0.03781513 -0.04201681 -0.04201681
[19,] 0.03781513 0.03781513 0.00000000 -0.08582488 0.00000000 0.00000000 0.00000000
[20,] 0.00000000 0.00000000 0.03781513 0.00000000 -0.08582488 0.00000000 0.00000000
[,15] [,16] [,17] [,18] [,19] [,20]
[1,] 0.00000000 -0.01019462 0.00000000 0.00000000 -0.01019462 -0.01019462
[2,] 0.00000000 0.00000000 -0.01019462 -0.01019462 0.00000000 0.00000000
[3,] 0.00000000 0.03707136 0.00000000 0.00000000 0.03707136 0.03707136
[4,] 0.00000000 0.00000000 0.03707136 0.03707136 0.00000000 0.00000000
[5,] 0.00000000 0.02521008 0.00000000 0.03781513 0.00000000 0.00000000
[6,] -0.01890756 0.00000000 0.03781513 0.03781513 0.00000000 0.00000000
[7,] -0.07563025 0.00000000 0.03781513 -0.01890756 0.00000000 0.03781513
[8,] 0.00000000 -0.07952236 -0.01890756 0.00000000 0.03781513 0.00000000
[9,] 0.03781513 -0.01890756 -0.10473244 0.00000000 0.03781513 0.00000000
[10,] -0.01890756 0.00000000 0.00000000 -0.10473244 0.00000000 0.03781513
[11,] 0.00000000 0.03781513 0.03781513 0.00000000 -0.08582488 0.00000000
[12,] 0.03781513 0.00000000 0.00000000 0.03781513 0.00000000 -0.08582488
[13,] 0.00000000 -0.02801120 0.00000000 -0.04201681 0.00000000 0.00000000
[14,] 0.02100840 0.00000000 -0.04201681 -0.04201681 0.00000000 0.00000000
[15,] 0.08403361 0.00000000 -0.04201681 0.02100840 0.00000000 -0.04201681
[16,] 0.00000000 0.11410217 0.02100840 0.00000000 -0.04201681 0.00000000
[17,] -0.04201681 0.02100840 0.14211338 0.00000000 -0.04201681 0.00000000
[18,] 0.02100840 0.00000000 0.00000000 0.14211338 0.00000000 -0.04201681
[19,] 0.00000000 -0.04201681 -0.04201681 0.00000000 0.12110498 0.00000000
[20,] -0.04201681 0.00000000 0.00000000 -0.04201681 0.00000000 0.12110498
>
> # RHS construction
> XPy = t(X) %*% Rif %*% y
> ZPy = t(Z) %*% Rif %*% y
>
> # print
> XPy
[,1]
[1,] 0.15449490
[2,] 0.06876738
[3,] 0.62001854
[4,] 0.36811863
> ZPy
[,1]
[1,] 0.00000000
[2,] 0.00000000
[3,] 0.00000000
[4,] 0.05579240
[5,] 0.02965709
[6,] 0.03911029
[7,] 0.03614458
[8,] 0.06255792
[9,] 0.00000000
[10,] 0.00000000
[11,] 0.00000000
[12,] 0.20620945
[13,] 0.15579240
[14,] 0.21232623
[15,] 0.18674699
[16,] 0.22706209
>
> RHS = rbind(XPy,
+ ZPy
+ )
>
> # print
> RHS
[,1]
[1,] 0.15449490
[2,] 0.06876738
[3,] 0.62001854
[4,] 0.36811863
[5,] 0.00000000
[6,] 0.00000000
[7,] 0.00000000
[8,] 0.05579240
[9,] 0.02965709
[10,] 0.03911029
[11,] 0.03614458
[12,] 0.06255792
[13,] 0.00000000
[14,] 0.00000000
[15,] 0.00000000
[16,] 0.20620945
[17,] 0.15579240
[18,] 0.21232623
[19,] 0.18674699
[20,] 0.22706209
>
> # Solutions
>
> # generalized inverse of LHS
> gi_LHS = ginv(LHS)
>
> # print
> gi_LHS
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 23.805131 6.291201 13.048276 5.671170 -6.5395332 -3.3374389 -5.2755138 -10.541038
[2,] 6.291201 31.996882 5.671170 16.041950 -5.3702695 -9.5573818 -4.3877316 -3.712060
[3,] 13.048276 5.671170 30.857577 12.600581 -5.8215059 -2.9950329 -4.8830435 -9.361379
[4,] 5.671170 16.041950 12.600581 38.483277 -4.9970359 -8.3727339 -3.6798298 -3.924119
[5,] -6.539533 -5.370269 -5.821506 -4.997036 18.6047347 0.3243012 1.5791386 8.540106
[6,] -3.337439 -9.557382 -2.995033 -8.372734 0.3243012 19.5963751 -0.4802746 1.001548
[7,] -5.275514 -4.387732 -4.883044 -3.679830 1.5791386 -0.4802746 17.8928862 2.301944
[8,] -10.541038 -3.712060 -9.361379 -3.924119 8.5401057 1.0015480 2.3019440 16.504656
[9,] -5.951896 -11.844286 -5.412983 -10.352380 2.1631118 9.3897226 7.6551424 2.269591
[10,] -6.630505 -12.149478 -5.929358 -10.731520 8.5774271 9.7250410 1.1203208 5.154530
[11,] -11.291251 -7.376935 -9.968148 -6.529986 5.6324645 4.7003325 5.0739444 9.706120
[12,] -9.583104 -7.784607 -8.815299 -6.559406 5.4460295 4.3104363 8.4506531 5.412337
[13,] -5.821506 -4.997036 -12.950918 -11.068126 16.0904195 0.4332042 2.1649074 7.012413
[14,] -2.995033 -8.372734 -6.657554 -18.656944 0.4332042 17.4241779 -0.6389804 1.388999
[15,] -4.883044 -3.679830 -10.821190 -8.237206 2.1649074 -0.6389804 15.1103333 3.128063
[16,] -9.361379 -3.924119 -20.830966 -8.590649 7.0124132 1.3889995 3.1280633 13.211993
[17,] -5.412983 -10.352380 -12.016345 -23.073618 2.9533095 8.1651708 5.7955284 3.105642
[18,] -5.929358 -10.731520 -13.184816 -23.892936 7.0407623 8.5802970 1.5641312 4.742595
[19,] -9.968148 -6.529986 -22.194546 -14.535359 5.3386160 4.0369934 4.6333538 8.655863
[20,] -8.815299 -6.559406 -19.547220 -14.675734 5.1134886 3.5591059 6.8877136 6.216282
[,9] [,10] [,11] [,12] [,13] [,14] [,15]
[1,] -5.951896 -6.630505 -11.291251 -9.583104 -5.8215059 -2.9950329 -4.8830435
[2,] -11.844286 -12.149478 -7.376935 -7.784607 -4.9970359 -8.3727339 -3.6798298
[3,] -5.412983 -5.929358 -9.968148 -8.815299 -12.9509185 -6.6575537 -10.8211898
[4,] -10.352380 -10.731520 -6.529986 -6.559406 -11.0681257 -18.6569441 -8.2372059
[5,] 2.163112 8.577427 5.632464 5.446029 16.0904195 0.4332042 2.1649074
[6,] 9.389723 9.725041 4.700333 4.310436 0.4332042 17.4241779 -0.6389804
[7,] 7.655142 1.120321 5.073944 8.450653 2.1649074 -0.6389804 15.1103333
[8,] 2.269591 5.154530 9.706120 5.412337 7.0124132 1.3889995 3.1280633
[9,] 16.541066 7.147506 8.548265 7.037833 2.9533095 8.1651708 5.7955284
[10,] 7.147506 17.151450 6.205604 8.531381 7.0407623 8.5802970 1.5641312
[11,] 8.548265 6.205604 17.115040 7.052592 5.3386160 4.0369934 4.6333538
[12,] 7.037833 8.531381 7.052592 16.284382 5.1134886 3.5591059 6.8877136
[13,] 2.953309 7.040762 5.338616 5.113489 35.9017859 0.9312686 4.6456425
[14,] 8.165171 8.580297 4.036993 3.559106 0.9312686 38.7676309 -1.3740157
[15,] 5.795528 1.564131 4.633354 6.887714 4.6456425 -1.3740157 33.7992438
[16,] 3.105642 4.742595 8.655863 6.216282 15.7328475 2.9783085 6.7165153
[17,] 13.278202 7.426558 7.024588 6.108719 6.3392413 18.2082861 13.1220853
[18,] 7.426558 14.036482 6.035384 7.010094 15.7970100 19.1056022 3.3523265
[19,] 7.024588 6.035384 13.970273 7.278308 11.8037249 9.0140533 10.2814964
[20,] 6.108719 7.010094 7.278308 12.951308 11.3161830 7.9802991 15.4655577
[,16] [,17] [,18] [,19] [,20]
[1,] -9.361379 -5.412983 -5.929358 -9.968148 -8.815299
[2,] -3.924119 -10.352380 -10.731520 -6.529986 -6.559406
[3,] -20.830966 -12.016345 -13.184816 -22.194546 -19.547220
[4,] -8.590649 -23.073618 -23.892936 -14.535359 -14.675734
[5,] 7.012413 2.953309 7.040762 5.338616 5.113489
[6,] 1.388999 8.165171 8.580297 4.036993 3.559106
[7,] 3.128063 5.795528 1.564131 4.633354 6.887714
[8,] 13.211993 3.105642 4.742595 8.655863 6.216282
[9,] 3.105642 13.278202 7.426558 7.024588 6.108719
[10,] 4.742595 7.426558 14.036482 6.035384 7.010094
[11,] 8.655863 7.024588 6.035384 13.970273 7.278308
[12,] 6.216282 6.108719 7.010094 7.278308 12.951308
[13,] 15.732847 6.339241 15.797010 11.803725 11.316183
[14,] 2.978309 18.208286 19.105602 9.014053 7.980299
[15,] 6.716515 13.122085 3.352327 10.281496 15.465558
[16,] 29.724917 6.665202 10.516097 19.252951 13.515032
[17,] 6.665202 29.864618 16.282618 15.758829 13.625005
[18,] 10.516097 16.282618 31.503255 13.311889 15.726464
[19,] 19.252951 15.758829 13.311889 31.363554 15.967135
[20,] 13.515032 13.625005 15.726464 15.967135 29.159492
>
> # solution
> sol = gi_LHS %*% RHS
>
> # print
> sol
[,1]
[1,] 4.360866999
[2,] 3.397261592
[3,] 6.799897620
[4,] 5.880295937
[5,] 0.150915567
[6,] -0.015392510
[7,] -0.078391896
[8,] -0.010238959
[9,] -0.270331441
[10,] 0.275808258
[11,] -0.316117562
[12,] 0.243755523
[13,] 0.279597973
[14,] -0.007610071
[15,] -0.170341439
[16,] -0.012670709
[17,] -0.477830262
[18,] 0.517238387
[19,] -0.478983695
[20,] 0.391961543
>
> # solutions for fixed effects and animal effects
> sol_t1_f1 = as.matrix(sol[1 : no_lvl_f1])
> sol_t2_f1 = as.matrix(sol[(no_lvl_f1 + 1) : (no_lvl_f1 * 2)])
> sol_f1 = cbind(sol_t1_f1, sol_t2_f1)
>
> #print
> sol_f1
[,1] [,2]
[1,] 4.360867 6.799898
[2,] 3.397262 5.880296
>
> # breedinv value
> sol_t1_bv = sol[(no_lvl_f1 * 2 + 1) : (no_lvl_f1 * 2 + no_animals)]
> sol_t2_bv = sol[(no_lvl_f1 * 2 + no_animals + 1) : (no_lvl_f1 * 2 + no_animals *2)]
> sol_bv = cbind(sol_t1_bv, sol_t2_bv)
>
> # print
> sol_bv
sol_t1_bv sol_t2_bv
[1,] 0.15091557 0.279597973
[2,] -0.01539251 -0.007610071
[3,] -0.07839190 -0.170341439
[4,] -0.01023896 -0.012670709
[5,] -0.27033144 -0.477830262
[6,] 0.27580826 0.517238387
[7,] -0.31611756 -0.478983695
[8,] 0.24375552 0.391961543
>
> # reliability(r2), accuracy(r), standard error of prediction(SEP)
>
> # diagonal elements of the generalized inverse of LHS for animal equation
> d_ani_t1 = diag(gi_LHS[(no_lvl_f1 * 2 + 1) : (no_lvl_f1 * 2 + no_animals),
+ (no_lvl_f1 * 2 + 1) : (no_lvl_f1 * 2 + no_animals)])
> d_ani_t2 = diag(gi_LHS[(no_lvl_f1 * 2 + no_animals + 1) : (no_lvl_f1 * 2 + no_animals * 2),
+ (no_lvl_f1 * 2 + no_animals + 1) : (no_lvl_f1 * 2 + no_animals * 2)])
> d_ani = cbind(d_ani_t1, d_ani_t2)
>
> # print
> d_ani
d_ani_t1 d_ani_t2
[1,] 18.60473 35.90179
[2,] 19.59638 38.76763
[3,] 17.89289 33.79924
[4,] 16.50466 29.72492
[5,] 16.54107 29.86462
[6,] 17.15145 31.50325
[7,] 17.11504 31.36355
[8,] 16.28438 29.15949
>
> # reliability
> rel = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
>
> for(i in 1 : no_animals){
+ rel[i, ] = 1 - d_ani[i,] / diag(G)
+ }
>
> # print
> rel
[,1] [,2]
[1,] 0.06976326 0.10245535
[2,] 0.02018124 0.03080923
[3,] 0.10535569 0.15501891
[4,] 0.17476719 0.25687708
[5,] 0.17294668 0.25338456
[6,] 0.14242750 0.21241863
[7,] 0.14424801 0.21591115
[8,] 0.18578089 0.27101269
>
> # accuracy
> acc = sqrt(rel)
>
> # print
> acc
[,1] [,2]
[1,] 0.2641274 0.3200865
[2,] 0.1420607 0.1755256
[3,] 0.3245854 0.3937244
[4,] 0.4180517 0.5068304
[5,] 0.4158686 0.5033732
[6,] 0.3773957 0.4608890
[7,] 0.3798000 0.4646624
[8,] 0.4310231 0.5205888
>
> # standard error of prediction(SEP)
> sep = sqrt(d_ani)
>
> # 확인
> sep
d_ani_t1 d_ani_t2
[1,] 4.313321 5.991810
[2,] 4.426779 6.226366
[3,] 4.229998 5.813712
[4,] 4.062592 5.452056
[5,] 4.067071 5.464853
[6,] 4.141431 5.612776
[7,] 4.137033 5.600317
[8,] 4.035391 5.399953
>
> # partitioning of breeding values
>
> # yield deviation
> yd1 = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
>
> # numerator of n2
> a2 = array(rep(0, no_of_trts * no_of_trts * no_animals), dim = c(no_of_trts, no_of_trts, no_animals))
>
> # looping data
> for (i in 1 : no_data) {
+ yd1[data[i, 1], ] = as.matrix(data[i, c(3,4)] - sol_f1[data[i,2],])
+
+ a2[,,data[i, 1]] = Ri
+ }
>
> # print
> yd1
[,1] [,2]
[1,] 0.0000000 0.0000000000
[2,] 0.0000000 0.0000000000
[3,] 0.0000000 0.0000000000
[4,] 0.1391330 0.0001023796
[5,] -0.4972616 -0.8802959373
[6,] 0.5027384 0.9197040627
[7,] -0.8608670 -0.7998976204
[8,] 0.6391330 0.7001023796
> a2
, , 1
[,1] [,2]
[1,] 0 0
[2,] 0 0
, , 2
[,1] [,2]
[1,] 0 0
[2,] 0 0
, , 3
[,1] [,2]
[1,] 0 0
[2,] 0 0
, , 4
[,1] [,2]
[1,] 0.02780352 -0.01019462
[2,] -0.01019462 0.03707136
, , 5
[,1] [,2]
[1,] 0.02780352 -0.01019462
[2,] -0.01019462 0.03707136
, , 6
[,1] [,2]
[1,] 0.02780352 -0.01019462
[2,] -0.01019462 0.03707136
, , 7
[,1] [,2]
[1,] 0.02780352 -0.01019462
[2,] -0.01019462 0.03707136
, , 8
[,1] [,2]
[1,] 0.02780352 -0.01019462
[2,] -0.01019462 0.03707136
>
> # Parents average, progeny contribution
>
> # parents avearge
> pa1 = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
>
> # progeny contribution numerator
> pc0 = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
>
> # numerator of n3, denominator of progeny contribution
> a3 = array(rep(0, no_of_trts * no_of_trts * no_animals), dim = c(no_of_trts, no_of_trts, no_animals))
>
> # numerator of n1
> a1 = array(rep(0, no_of_trts * no_of_trts * no_animals), dim = c(no_of_trts, no_of_trts, no_animals))
>
>
> # looping pedi
> for (i in 1 : no_animals) {
+
+ sire = pedi[i, 2]
+ dam = pedi[i, 3]
+
+ if (sire == 0 & dam == 0) {
+ # both parents unknown
+ # PA
+ a1[ , , i] = 1 * Gi
+
+ } else if (sire != 0 & dam == 0) {
+ # 개체의 부만 알고 있을 경우
+
+ # PA
+ a1[ , , i] = 4/3 * Gi
+ pa1[i, ] = sol_bv[sire, ]/2
+
+ # PC for sire
+ a3[ , , sire] = a3[ , , sire] + 0.5 * Gi * (2/3)
+ pc0[sire, ] = pc0[sire, ] + (0.5 * Gi * (2/3)) %*% (2 * sol_bv[i,])
+
+ } else if (sire == 0 & dam != 0) {
+ # dam known
+
+ # PA
+ a1[ , , i] = 4/3 * Gi
+ pa1[i, ] = sol_bv[dam, ]/2
+
+ # PC for dam
+ a3[ , , dam] = a3[ , , dam] + 0.5 * Gi * (2/3)
+ pc0[dam, ] = pc0[dam, ] + (0.5 * Gi * (2/3)) %*% (2 * sol_bv[i])
+
+ } else {
+ # both parents known
+
+ # PA
+ a1[ , , i] = 2 * Gi
+ pa1[i, ] = (sol_bv[sire, ] + sol_bv[dam,])/2
+
+ # PC for sire
+ a3[ , , sire] = a3[ , , sire] + 0.5 * Gi
+ pc0[sire, ] = pc0[sire, ] + (0.5 * Gi) %*% (2 * sol_bv[i, ] - sol_bv[dam, ])
+
+ # PC for dam
+ a3[ , , dam] = a3[ , , dam] + 0.5 * Gi
+ pc0[dam, ] = pc0[dam, ] + (0.5 * Gi) %*% (2 * sol_bv[i, ] - sol_bv[sire, ])
+
+ }
+ }
>
> # print
> a1
, , 1
[,1] [,2]
[1,] 0.08403361 -0.03781513
[2,] -0.03781513 0.04201681
, , 2
[,1] [,2]
[1,] 0.08403361 -0.03781513
[2,] -0.03781513 0.04201681
, , 3
[,1] [,2]
[1,] 0.08403361 -0.03781513
[2,] -0.03781513 0.04201681
, , 4
[,1] [,2]
[1,] 0.11204482 -0.05042017
[2,] -0.05042017 0.05602241
, , 5
[,1] [,2]
[1,] 0.16806723 -0.07563025
[2,] -0.07563025 0.08403361
, , 6
[,1] [,2]
[1,] 0.16806723 -0.07563025
[2,] -0.07563025 0.08403361
, , 7
[,1] [,2]
[1,] 0.16806723 -0.07563025
[2,] -0.07563025 0.08403361
, , 8
[,1] [,2]
[1,] 0.16806723 -0.07563025
[2,] -0.07563025 0.08403361
> pa1
[,1] [,2]
[1,] 0.00000000 0.00000000
[2,] 0.00000000 0.00000000
[3,] 0.00000000 0.00000000
[4,] 0.07545778 0.13979899
[5,] -0.04689220 -0.08897575
[6,] 0.06776153 0.13599395
[7,] -0.14028520 -0.24525049
[8,] 0.09870818 0.17344847
> a3
, , 1
[,1] [,2]
[1,] 0.07002801 -0.03151261
[2,] -0.03151261 0.03501401
, , 2
[,1] [,2]
[1,] 0.08403361 -0.03781513
[2,] -0.03781513 0.04201681
, , 3
[,1] [,2]
[1,] 0.08403361 -0.03781513
[2,] -0.03781513 0.04201681
, , 4
[,1] [,2]
[1,] 0.04201681 -0.01890756
[2,] -0.01890756 0.02100840
, , 5
[,1] [,2]
[1,] 0.04201681 -0.01890756
[2,] -0.01890756 0.02100840
, , 6
[,1] [,2]
[1,] 0.04201681 -0.01890756
[2,] -0.01890756 0.02100840
, , 7
[,1] [,2]
[1,] 0 0
[2,] 0 0
, , 8
[,1] [,2]
[1,] 0 0
[2,] 0 0
> pc0
[,1] [,2]
[1,] 0.0038664044 0.0110750251
[2,] -0.0020114248 0.0005246376
[3,] -0.0002921427 -0.0083856077
[4,] -0.0061278140 -0.0032441978
[5,] -0.0082610361 -0.0080987423
[6,] 0.0057346179 0.0093477285
[7,] 0.0000000000 0.0000000000
[8,] 0.0000000000 0.0000000000
>
> # denominator of n1, n2, n3, diagonal of animals in LHS
> n_de = a1 + a2 + a3
>
> # print
> n_de
, , 1
[,1] [,2]
[1,] 0.15406162 -0.06932773
[2,] -0.06932773 0.07703081
, , 2
[,1] [,2]
[1,] 0.16806723 -0.07563025
[2,] -0.07563025 0.08403361
, , 3
[,1] [,2]
[1,] 0.16806723 -0.07563025
[2,] -0.07563025 0.08403361
, , 4
[,1] [,2]
[1,] 0.18186515 -0.07952236
[2,] -0.07952236 0.11410217
, , 5
[,1] [,2]
[1,] 0.2378876 -0.1047324
[2,] -0.1047324 0.1421134
, , 6
[,1] [,2]
[1,] 0.2378876 -0.1047324
[2,] -0.1047324 0.1421134
, , 7
[,1] [,2]
[1,] 0.19587075 -0.08582488
[2,] -0.08582488 0.12110498
, , 8
[,1] [,2]
[1,] 0.19587075 -0.08582488
[2,] -0.08582488 0.12110498
>
> # parents average fraction of breeding values
> pa = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
> for (i in 1 : no_animals) {
+ pa[i, ] = ginv(n_de[ , , i]) %*% a1[, , i] %*% pa1[i, ]
+ }
>
> # print
> pa
[,1] [,2]
[1,] 0.00000000 0.00000000
[2,] 0.00000000 0.00000000
[3,] 0.00000000 0.00000000
[4,] 0.03331733 0.05851558
[5,] -0.02519179 -0.04622283
[6,] 0.03577082 0.07071542
[7,] -0.08971192 -0.14614589
[8,] 0.06301831 0.10337078
>
> # yield deviation fraction of breeding values
> yd = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
> for (i in 1 : no_animals) {
+ yd[i, ] = ginv(n_de[ , , i]) %*% a2[, , i] %*% yd1[i, ]
+ }
>
> # print
> yd
[,1] [,2]
[1,] 0.0000000 0.000000000
[2,] 0.0000000 0.000000000
[3,] 0.0000000 0.000000000
[4,] 0.0227885 0.003484439
[5,] -0.1565945 -0.309364950
[6,] 0.1614856 0.322856538
[7,] -0.2264056 -0.332837809
[8,] 0.1807372 0.288590763
>
> # progeny contribution
> pc1 = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
> for (i in 1 : no_animals) {
+ pc1[i, ] = ginv(a3[ , , i]) %*% pc0[i, ]
+ }
> pc1[is.nan(pc1) == TRUE] = 0
> pc1
[,1] [,2]
[1,] 0.33201425 0.61511554
[2,] -0.03078502 -0.01522014
[3,] -0.15678379 -0.34068288
[4,] -0.36190368 -0.48013713
[5,] -0.62199616 -0.94529668
[6,] 0.56590294 0.95426452
[7,] 0.00000000 0.00000000
[8,] 0.00000000 0.00000000
>
> # progeny contribution fraction of breeding values
> pc = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
> for (i in 1 : no_animals) {
+ pc[i, ] = ginv(n_de[ , , i]) %*% a3[, , i] %*% pc1[i, ]
+ }
>
> # print
> pc
[,1] [,2]
[1,] 0.15091557 0.279597973
[2,] -0.01539251 -0.007610071
[3,] -0.07839190 -0.170341439
[4,] -0.06634480 -0.074670727
[5,] -0.08854515 -0.122242479
[6,] 0.07855184 0.123666429
[7,] 0.00000000 0.000000000
[8,] 0.00000000 0.000000000
>
> # Progeny(Daughter) Yield Deviation(PYD, DYD)
>
> # n2 of progeny
> n2prog = array(rep(0, no_of_trts * no_of_trts * no_animals), dim = c(no_of_trts, no_of_trts, no_animals))
> for (i in 1 : no_animals) {
+ n2prog[ , , i] = ginv((a1 + a2)[, , i]) %*% a2[, , i]
+ }
>
> # print
> n2prog
, , 1
[,1] [,2]
[1,] 0 0
[2,] 0 0
, , 2
[,1] [,2]
[1,] 0 0
[2,] 0 0
, , 3
[,1] [,2]
[1,] 0 0
[2,] 0 0
, , 4
[,1] [,2]
[1,] 0.21085286 0.1389018
[2,] 0.02778035 0.4886564
, , 5
[,1] [,2]
[1,] 0.15238095 0.1190476
[2,] 0.02380952 0.3904762
, , 6
[,1] [,2]
[1,] 0.15238095 0.1190476
[2,] 0.02380952 0.3904762
, , 7
[,1] [,2]
[1,] 0.15238095 0.1190476
[2,] 0.02380952 0.3904762
, , 8
[,1] [,2]
[1,] 0.15238095 0.1190476
[2,] 0.02380952 0.3904762
>
> # numerator of dyd : summation of u of progeny * n2 of progeny * (2 * YD - bv of mate)
> dyd_n = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
> # denominator of dyd : summation of u of progeny * n2 of progeny
> dyd_d = array(rep(0, no_of_trts * no_of_trts * no_animals), dim = c(no_of_trts, no_of_trts, no_animals))
>
> # looping pedi
> for (i in 1 : no_animals) {
+ # i = 5
+ sire = pedi[i, 2]
+ dam = pedi[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_bv[dam, ])
+ dyd_n[dam, ] = dyd_n[dam, ] + n2prog[ , , i] %*% (2 * yd1[i, ] - sol_bv[sire, ])
+
+ # dyd_d
+ dyd_d[ , , sire] = dyd_d[ , , sire] + n2prog[ , , i]
+ dyd_d[ , , dam] = dyd_d[ , , dam] + n2prog[ , , i]
+
+ }
+ }
>
> # print
> dyd_n
[,1] [,2]
[1,] 0.41457857 0.75074330
[2,] -0.01300594 -0.01335216
[3,] -0.10001866 -0.33916490
[4,] -0.35473336 -0.47265781
[5,] -0.44974264 -0.66048422
[6,] 0.39369861 0.64556227
[7,] 0.00000000 0.00000000
[8,] 0.00000000 0.00000000
> dyd_d
, , 1
[,1] [,2]
[1,] 0.29294952 0.2116488
[2,] 0.04232976 0.7162471
, , 2
[,1] [,2]
[1,] 0.30476190 0.2380952
[2,] 0.04761905 0.7809524
, , 3
[,1] [,2]
[1,] 0.30476190 0.2380952
[2,] 0.04761905 0.7809524
, , 4
[,1] [,2]
[1,] 0.15238095 0.1190476
[2,] 0.02380952 0.3904762
, , 5
[,1] [,2]
[1,] 0.15238095 0.1190476
[2,] 0.02380952 0.3904762
, , 6
[,1] [,2]
[1,] 0.15238095 0.1190476
[2,] 0.02380952 0.3904762
, , 7
[,1] [,2]
[1,] 0 0
[2,] 0 0
, , 8
[,1] [,2]
[1,] 0 0
[2,] 0 0
>
> # dyd
> dyd = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
>
> # looping pedi
> for (i in 1 : no_animals) {
+ dyd[i, ] = ginv(dyd_d[ , , i]) %*% dyd_n[i, ]
+ }
> dyd[is.nan(dyd) == TRUE] = 0
>
> # print
> dyd
[,1] [,2]
[1,] 0.68726086 1.00754574
[2,] -0.03078502 -0.01522014
[3,] 0.01166354 -0.43500772
[4,] -1.45140256 -1.12196498
[5,] -1.71149504 -1.58712453
[6,] 1.35665790 1.57054620
[7,] 0.00000000 0.00000000
[8,] 0.00000000 0.00000000
>
> # breeding values and fractions
> mt1_result = data.frame(animal = pedi[,1], animal_bv = sol_bv, rel = rel, acc = acc, sep = sep, pa = pa, yd = yd, pc = pc, sum_of_fr = pa + yd + pc, dyd = dyd)
>
> # print
> mt1_result
animal animal_bv.sol_t1_bv animal_bv.sol_t2_bv rel.1 rel.2 acc.1 acc.2
1 1 0.15091557 0.279597973 0.06976326 0.10245535 0.2641274 0.3200865
2 2 -0.01539251 -0.007610071 0.02018124 0.03080923 0.1420607 0.1755256
3 3 -0.07839190 -0.170341439 0.10535569 0.15501891 0.3245854 0.3937244
4 4 -0.01023896 -0.012670709 0.17476719 0.25687708 0.4180517 0.5068304
5 5 -0.27033144 -0.477830262 0.17294668 0.25338456 0.4158686 0.5033732
6 6 0.27580826 0.517238387 0.14242750 0.21241863 0.3773957 0.4608890
7 7 -0.31611756 -0.478983695 0.14424801 0.21591115 0.3798000 0.4646624
8 8 0.24375552 0.391961543 0.18578089 0.27101269 0.4310231 0.5205888
sep.d_ani_t1 sep.d_ani_t2 pa.1 pa.2 yd.1 yd.2 pc.1
1 4.313321 5.991810 0.00000000 0.00000000 0.0000000 0.000000000 0.15091557
2 4.426779 6.226366 0.00000000 0.00000000 0.0000000 0.000000000 -0.01539251
3 4.229998 5.813712 0.00000000 0.00000000 0.0000000 0.000000000 -0.07839190
4 4.062592 5.452056 0.03331733 0.05851558 0.0227885 0.003484439 -0.06634480
5 4.067071 5.464853 -0.02519179 -0.04622283 -0.1565945 -0.309364950 -0.08854515
6 4.141431 5.612776 0.03577082 0.07071542 0.1614856 0.322856538 0.07855184
7 4.137033 5.600317 -0.08971192 -0.14614589 -0.2264056 -0.332837809 0.00000000
8 4.035391 5.399953 0.06301831 0.10337078 0.1807372 0.288590763 0.00000000
pc.2 sum_of_fr.1 sum_of_fr.2 dyd.1 dyd.2
1 0.279597973 0.15091557 0.279597973 0.68726086 1.00754574
2 -0.007610071 -0.01539251 -0.007610071 -0.03078502 -0.01522014
3 -0.170341439 -0.07839190 -0.170341439 0.01166354 -0.43500772
4 -0.074670727 -0.01023896 -0.012670709 -1.45140256 -1.12196498
5 -0.122242479 -0.27033144 -0.477830262 -1.71149504 -1.58712453
6 0.123666429 0.27580826 0.517238387 1.35665790 1.57054620
7 0.000000000 -0.31611756 -0.478983695 0.00000000 0.00000000
8 0.000000000 0.24375552 0.391961543 0.00000000 0.00000000
>
> # 파일로 출력, 분리자는 ",", 따옴표 없고
> output_filename = gsub("[ -]", "", paste("mt1_result_", Sys.Date(), ".csv"))
> write.table(mt1_result, output_filename, sep=",", row.names = FALSE, quote = FALSE)
>
다음은 출력한 파일이다.
animal,animal_bv.sol_t1_bv,animal_bv.sol_t2_bv,rel.1,rel.2,acc.1,acc.2,sep.d_ani_t1,sep.d_ani_t2,pa.1,pa.2,yd.1,yd.2,pc.1,pc.2,sum_of_fr.1,sum_of_fr.2,dyd.1,dyd.2
1,0.150915567312797,0.279597972744865,0.0697632637735797,0.102455353618367,0.264127362788447,0.320086478343536,4.31332061462261,5.99180989812472,0,0,0,0,0.150915567312795,0.279597972744862,0.150915567312795,0.279597972744862,0.68726085711682,1.00754574375344
2,-0.0153925096914798,-0.0076100708199065,0.0201812427467539,0.030809227182284,0.14206070092307,0.175525574154549,4.42677931967078,6.22636578693451,0,0,0,0,-0.0153925096914804,-0.00761007081990604,-0.0153925096914804,-0.00761007081990604,-0.0307850193829532,-0.015220141639821
3,-0.0783918961641927,-0.170341438593948,0.105355689244771,0.155018905719144,0.324585411324617,0.39372440325581,4.22999837057942,5.81371170348464,0,0,0,0,-0.078391896164192,-0.170341438593946,-0.078391896164192,-0.170341438593946,0.0116635350966923,-0.435007715776464
4,-0.0102389585292884,-0.0126707086241291,0.174767187326832,0.256877083095507,0.418051656290024,0.506830428344143,4.06259230706988,5.45205618791477,0.0333173344095044,0.0585155786066336,0.0227885035278452,0.00348443941480283,-0.0663447964666395,-0.0746707266455604,-0.0102389585292899,-0.0126707086241239,-1.45140255671963,-1.12196497916736
5,-0.270331441389235,-0.477830261602733,0.172946680277051,0.253384555790556,0.415868585345239,0.50337317746435,4.06707098468406,5.46485295029773,-0.0251917915640488,-0.0462228319310023,-0.156594500241952,-0.3093649501734,-0.0885451495832349,-0.122242479498329,-0.270331441389235,-0.477830261602731,-1.71149503957958,-1.58712453214596
6,0.275808257580577,0.517238387038379,0.142427497991221,0.212418627668353,0.377395678289009,0.460888953727851,4.14143091698698,5.61277604160952,0.0357708240803181,0.070715419872911,0.161485595241754,0.322856537923106,0.078551838258502,0.123666429242362,0.275808257580574,0.517238387038379,1.35665789805533,1.57054619782385
7,-0.316117561639437,-0.478983695055088,0.144248005041,0.215911154973304,0.379799953977091,0.46466240968396,4.13703274088809,5.60031729467785,-0.0897119212614896,-0.146145886165347,-0.226405640377943,-0.332837808889743,0,0,-0.316117561639433,-0.47898369505509,0,0
8,0.2437555230054,0.391961542524077,0.185780891703132,0.271012693707549,0.431023075604001,0.52058879521898,4.03539120358081,5.39995298606368,0.0630183062404894,0.103370779985251,0.180737216764914,0.288590762538829,0,0,0.243755523005404,0.391961542524079,0,0
엑셀에서 열면 다음과 같다.
