반복 모형은 한 개체에 대하여 긴 시간 동안 여러번 능력을 조사(측정)한 것을 분석하는 모형이다. 예를 들어 육우에 대해서 6개월령, 12개월령, 18개월령 체중을 잴 수 있다. 젖소의 경우 1산차, 2산차 유량을 생산할 수 있다. 이런 자료를 반복 모형으로 분석할 때의 이론적 배경은 생략한다.
자료는 다음과 같다.
4 1 1 201
4 2 3 280
5 1 1 150
5 2 4 200
6 1 2 160
6 2 3 190
7 1 1 180
7 2 3 250
8 1 2 285
8 2 4 300
animal parity hys(herd-year-season) fat-yield이다.
repeat-data.txt로 저장
혈통은 다음과 같다.
1 0 0
2 0 0
3 0 0
4 1 2
5 3 2
6 1 5
7 3 4
8 1 7
animal sire dam이다.
pedi.txt로 저장한다.
육종가를 구하고, 신뢰도, 정확도, SEP를 구하고, 육종가를 분할하고, DYD를 구하는 R 프로그램은 다음과 같다.
# Repeatability Model - Date : 2020-06-14
# Linear Models for the Prediction of Animal Breeding Values, 3rd Edtion
# Raphael Mrode
# Example 4.1 p62
# install.packages('MASS', dependencies = TRUE)
library(MASS)
# 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/repeatability_r")
# print working directory
getwd()
# prdigree and data
# read pedigree : animal sire dam
pedi = read.table("pedi.txt", header = FALSE, sep = "", col.names = c("animal", "sire", "dam"),stringsAsFactors = FALSE)
# print
pedi
# read data : animal, parity, hys, fat_yield
data = read.table("repeat_data.txt", header = FALSE, sep = "", col.names = c("animal", "parity", "hys", "hatyield"), stringsAsFactors = FALSE)
# print
data
# variance component and ratio
sigma_a = 20
sigma_e = 28
sigma_pe = 12
alpha = sigma_e/sigma_a
alpha2 = sigma_e/sigma_pe
# print
sigma_a
sigma_e
sigma_pe
alpha
alpha2
# design matrix
# design matrix of hys fixed effect
X1 = desgn(data[, 3])
# print
X1
# number of levels of hys fixed levels
no_lvl_f1 = ncol(X1)
# print
no_lvl_f1
# design matrix of parity fixed effect
X2 = desgn(data[, 2])
# print
X2
# number of levels of parity fixed levels
no_lvl_f2 = ncol(X2)
# print
no_lvl_f2
# final fixed design matrix : X - hys, parity
X = cbind(X1, X2)
# print
X
# number of levels of fixed effect
no_lvl_f = ncol(X)
# print
no_lvl_f
# desing matrix of animal effect
Z = desgn(data[, 1])
# print
Z
# number of animals
no_animals = ncol(Z)
# print
no_animals
# observation
y = data[, 4]
# print
y
# inverse matrix of NRM
ai = ainv(pedi)
# print
ai
# LHS, RHS
# LHS construction
XPX = t(X) %*% X
XPZ = t(X) %*% Z
ZPX = t(Z) %*% X
ZPZ = t(Z) %*% Z
#print
XPX
XPZ
ZPX
ZPZ
LHS = rbind(
cbind(XPX, XPZ, XPZ),
cbind(ZPX, ZPZ + ai * alpha, ZPZ),
cbind(ZPX, ZPZ, ZPZ + diag(no_animals) * alpha2)
)
# print
LHS
# RHS construction
XPy = t(X) %*% y
ZPy = t(Z) %*% y
# print
XPy
RHS = rbind(
XPy,
ZPy,
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
# hys
sol_hys = sol[1 : no_lvl_f1]
# parity
sol_parity = sol[(no_lvl_f1 + 1) : (no_lvl_f1 + no_lvl_f2)]
# animal
sol_animal = sol[(no_lvl_f1 + no_lvl_f2 + 1) : (no_lvl_f1 + no_lvl_f2 + no_animals)]
# permanent environmental
sol_pe = sol[(no_lvl_f1 + no_lvl_f2 + no_animals + 1) : (no_lvl_f1 + no_lvl_f2 + no_animals * 2)]
# print
sol_hys
sol_parity
sol_animal
sol_pe
# reliability(r2), accuracy(r), standard error of prediction(SEP)
# diagonal elements of the generalized inverse of LHS for animal equation
d_ani = diag(gi_LHS[(no_lvl_f1 + no_lvl_f2 + 1) : (no_lvl_f1 + no_lvl_f2 + no_animals),
(no_lvl_f1 + no_lvl_f2 + 1) : (no_lvl_f1 + no_lvl_f2 + no_animals)])
# print
d_ani
# reliability
rel = 1 - d_ani * alpha
# print
rel
# accuracy
acc = sqrt(rel)
# print
acc
# standard error of prediction(SEP)
sep = sqrt(d_ani * sigma_e)
# 확인
sep
# partitioning of breeding values
# yield deviation
yd1 = ginv(ZPZ) %*% t(Z) %*% (y - X1 %*% sol_hys - X2 %*% sol_parity - Z %*% sol_pe)
# print
yd1
# numerator of n2
a2 = diag(ZPZ)
# print
a2
# Parents average, progeny contribution
# parents avearge
pa1 = rep(0, no_animals)
# progeny contribution numerator
pc0 = rep(0, no_animals)
# numerator of n3, denominator of progeny contribution
a3 = rep(0, no_animals)
# numerator of n1
a1 = rep(0, 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 * alpha
} else if (sire != 0 & dam == 0) {
# 개체의 부만 알고 있을 경우
# PA
a1[i] = 4/3 * alpha
pa1[i] = sol_animal[sire]/2
# PC for sire
a3[sire] = a3[sire] + 0.5 * alpha * (2/3)
pc0[sire] = pc0[sire] + 0.5 * alpha * (2/3) * (2 * sol_animal[i])
} else if (sire == 0 & dam != 0) {
# dam known
# PA
a1[i] = 4/3 * alpha
pa1[i] = sol_animal[dam]/2
# PC for dam
a3[dam] = a3[dam] + 0.5 * alpha * (2/3)
pc0[dam] = pc0[dam] + 0.5 * alpha * (2/3) * (2 * sol_animal[i])
} else {
# both parents known
# PA
a1[i] = 2 * alpha
pa1[i] = (sol_animal[sire] + sol_animal[dam])/2
# PC for sire
a3[sire] = a3[sire] + 0.5 * alpha
pc0[sire] = pc0[sire] + 0.5 * alpha * (2 * sol_animal[i] - sol_animal[dam])
# PC for dam
a3[dam] = a3[dam] + 0.5 * alpha
pc0[dam] = pc0[dam] + 0.5 * alpha * (2 * sol_animal[i] - sol_animal[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 = pa1 * (a1 / n_de)
# print
pa
# yield deviation fraction of breeding values
yd = yd1 * (a2 / n_de)
# print
yd
# progeny contribution
pc1 = pc0 / a3
pc1[is.nan(pc1) == TRUE] = 0
pc1
# progeny contribution fraction of breeding values
pc = pc1 * (a3 / n_de)
# print
pc
# Progeny(Daughter) Yield Deviation(PYD, DYD)
# n2 of progeny
n2prog = a2 / (a1 + a2)
# print
n2prog
# numerator of dyd : summation of u of progeny * n2 of progeny * (2 * YD - bv of mate)
dyd_n = rep(0, no_animals)
# denominator of dyd : summation of u of progeny * n2 of progeny
dyd_d = rep(0, 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
} 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]
}
}
# print
dyd_n
dyd_d
# dyd
dyd = dyd_n / dyd_d
dyd[is.nan(dyd) == TRUE] = 0
# print
dyd
# breeding values and fractions
repeat_result = data.frame(animal = pedi[,1], animal_bv = sol_animal, rel = rel, acc = acc, sep = sep, pa = pa, yd = yd, pc = pc, sum_of_fr = pa + yd + pc, dyd = dyd)
# print
repeat_result
# 파일로 출력, 분리자는 ",", 따옴표 없고
output_filename = gsub("[ -]", "", paste("repeat_result_", Sys.Date(), ".csv"))
write.table(repeat_result, output_filename, sep=",", row.names = FALSE, quote = FALSE)
다음은 실행 결과다.
# Repeatability Model - Date : 2020-06-14
>
> # Linear Models for the Prediction of Animal Breeding Values, 3rd Edtion
>
> # Raphael Mrode
>
> # Example 4.1 p62
>
> # install.packages('MASS', dependencies = TRUE)
> library(MASS)
>
> # 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/repeatability_r")
>
> # print working directory
> getwd()
[1] "D:/temp/repeatability_r"
>
>
> # prdigree and data
>
> # read pedigree : animal sire dam
> pedi = read.table("pedi.txt", header = FALSE, sep = "", col.names = c("animal", "sire", "dam"),stringsAsFactors = FALSE)
>
> # print
> pedi
animal sire dam
1 1 0 0
2 2 0 0
3 3 0 0
4 4 1 2
5 5 3 2
6 6 1 5
7 7 3 4
8 8 1 7
>
> # read data : animal, parity, hys, fat_yield
> data = read.table("repeat_data.txt", header = FALSE, sep = "", col.names = c("animal", "parity", "hys", "hatyield"), stringsAsFactors = FALSE)
>
> # print
> data
animal parity hys hatyield
1 4 1 1 201
2 4 2 3 280
3 5 1 1 150
4 5 2 4 200
5 6 1 2 160
6 6 2 3 190
7 7 1 1 180
8 7 2 3 250
9 8 1 2 285
10 8 2 4 300
>
> # variance component and ratio
> sigma_a = 20
> sigma_e = 28
> sigma_pe = 12
>
> alpha = sigma_e/sigma_a
> alpha2 = sigma_e/sigma_pe
>
> # print
> sigma_a
[1] 20
> sigma_e
[1] 28
> sigma_pe
[1] 12
> alpha
[1] 1.4
> alpha2
[1] 2.333333
>
>
> # design matrix
>
> # design matrix of hys fixed effect
> X1 = desgn(data[, 3])
>
> # print
> 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
>
> # number of levels of hys fixed levels
> no_lvl_f1 = ncol(X1)
>
> # print
> no_lvl_f1
[1] 4
>
>
> # design matrix of parity fixed effect
> X2 = desgn(data[, 2])
>
> # print
> 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
>
> # number of levels of parity fixed levels
> no_lvl_f2 = ncol(X2)
>
> # print
> no_lvl_f2
[1] 2
>
> # final fixed design matrix : X - hys, parity
> X = cbind(X1, X2)
>
> # print
> 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
>
> # number of levels of fixed effect
> no_lvl_f = ncol(X)
>
> # print
> no_lvl_f
[1] 6
>
> # desing matrix of animal effect
> Z = desgn(data[, 1])
>
> # print
> 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
>
> # number of animals
> no_animals = ncol(Z)
>
> # print
> no_animals
[1] 8
>
>
> # observation
> y = data[, 4]
>
> # print
> y
[1] 201 280 150 200 160 190 180 250 285 300
>
> # inverse matrix of NRM
> ai = ainv(pedi)
>
> # print
> ai
[,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
> XPX = t(X) %*% X
> XPZ = t(X) %*% Z
> ZPX = t(Z) %*% X
> ZPZ = t(Z) %*% Z
>
> #print
> XPX
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3 0 0 0 3 0
[2,] 0 2 0 0 2 0
[3,] 0 0 3 0 0 3
[4,] 0 0 0 2 0 2
[5,] 3 2 0 0 5 0
[6,] 0 0 3 2 0 5
> XPZ
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 0 0 0 1 1 0 1 0
[2,] 0 0 0 0 0 1 0 1
[3,] 0 0 0 1 0 1 1 0
[4,] 0 0 0 0 1 0 0 1
[5,] 0 0 0 1 1 1 1 1
[6,] 0 0 0 1 1 1 1 1
> ZPX
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0 0 0 0 0 0
[2,] 0 0 0 0 0 0
[3,] 0 0 0 0 0 0
[4,] 1 0 1 0 1 1
[5,] 1 0 0 1 1 1
[6,] 0 1 1 0 1 1
[7,] 1 0 1 0 1 1
[8,] 0 1 0 1 1 1
> ZPZ
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 0 0 0 0 0 0 0 0
[2,] 0 0 0 0 0 0 0 0
[3,] 0 0 0 0 0 0 0 0
[4,] 0 0 0 2 0 0 0 0
[5,] 0 0 0 0 2 0 0 0
[6,] 0 0 0 0 0 2 0 0
[7,] 0 0 0 0 0 0 2 0
[8,] 0 0 0 0 0 0 0 2
>
> LHS = rbind(
+ cbind(XPX, XPZ, XPZ),
+ cbind(ZPX, ZPZ + ai * alpha, ZPZ),
+ cbind(ZPX, ZPZ, ZPZ + diag(no_animals) * alpha2)
+ )
>
> # print
> LHS
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
[1,] 3 0 0 0 3 0 0.0 0.0 0.0 1.0 1.0 0.0 1.0 0.0 0.000000
[2,] 0 2 0 0 2 0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 1.0 0.000000
[3,] 0 0 3 0 0 3 0.0 0.0 0.0 1.0 0.0 1.0 1.0 0.0 0.000000
[4,] 0 0 0 2 0 2 0.0 0.0 0.0 0.0 1.0 0.0 0.0 1.0 0.000000
[5,] 3 2 0 0 5 0 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 0.000000
[6,] 0 0 3 2 0 5 0.0 0.0 0.0 1.0 1.0 1.0 1.0 1.0 0.000000
[7,] 0 0 0 0 0 0 3.5 0.7 0.0 -1.4 0.7 -1.4 0.7 -1.4 0.000000
[8,] 0 0 0 0 0 0 0.7 2.8 0.7 -1.4 -1.4 0.0 0.0 0.0 0.000000
[9,] 0 0 0 0 0 0 0.0 0.7 2.8 0.7 -1.4 0.0 -1.4 0.0 0.000000
[10,] 1 0 1 0 1 1 -1.4 -1.4 0.7 5.5 0.0 0.0 -1.4 0.0 0.000000
[11,] 1 0 0 1 1 1 0.7 -1.4 -1.4 0.0 5.5 -1.4 0.0 0.0 0.000000
[12,] 0 1 1 0 1 1 -1.4 0.0 0.0 0.0 -1.4 4.8 0.0 0.0 0.000000
[13,] 1 0 1 0 1 1 0.7 0.0 -1.4 -1.4 0.0 0.0 5.5 -1.4 0.000000
[14,] 0 1 0 1 1 1 -1.4 0.0 0.0 0.0 0.0 0.0 -1.4 4.8 0.000000
[15,] 0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.333333
[16,] 0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.000000
[17,] 0 0 0 0 0 0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.000000
[18,] 1 0 1 0 1 1 0.0 0.0 0.0 2.0 0.0 0.0 0.0 0.0 0.000000
[19,] 1 0 0 1 1 1 0.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 0.000000
[20,] 0 1 1 0 1 1 0.0 0.0 0.0 0.0 0.0 2.0 0.0 0.0 0.000000
[21,] 1 0 1 0 1 1 0.0 0.0 0.0 0.0 0.0 0.0 2.0 0.0 0.000000
[22,] 0 1 0 1 1 1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.0 0.000000
[,16] [,17] [,18] [,19] [,20] [,21] [,22]
[1,] 0.000000 0.000000 1.000000 1.000000 0.000000 1.000000 0.000000
[2,] 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000 1.000000
[3,] 0.000000 0.000000 1.000000 0.000000 1.000000 1.000000 0.000000
[4,] 0.000000 0.000000 0.000000 1.000000 0.000000 0.000000 1.000000
[5,] 0.000000 0.000000 1.000000 1.000000 1.000000 1.000000 1.000000
[6,] 0.000000 0.000000 1.000000 1.000000 1.000000 1.000000 1.000000
[7,] 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
[8,] 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
[9,] 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
[10,] 0.000000 0.000000 2.000000 0.000000 0.000000 0.000000 0.000000
[11,] 0.000000 0.000000 0.000000 2.000000 0.000000 0.000000 0.000000
[12,] 0.000000 0.000000 0.000000 0.000000 2.000000 0.000000 0.000000
[13,] 0.000000 0.000000 0.000000 0.000000 0.000000 2.000000 0.000000
[14,] 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.000000
[15,] 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
[16,] 2.333333 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
[17,] 0.000000 2.333333 0.000000 0.000000 0.000000 0.000000 0.000000
[18,] 0.000000 0.000000 4.333333 0.000000 0.000000 0.000000 0.000000
[19,] 0.000000 0.000000 0.000000 4.333333 0.000000 0.000000 0.000000
[20,] 0.000000 0.000000 0.000000 0.000000 4.333333 0.000000 0.000000
[21,] 0.000000 0.000000 0.000000 0.000000 0.000000 4.333333 0.000000
[22,] 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 4.333333
>
> # RHS construction
> XPy = t(X) %*% y
> ZPy = t(Z) %*% y
>
> # print
> XPy
[,1]
[1,] 531
[2,] 445
[3,] 720
[4,] 500
[5,] 976
[6,] 1220
>
> RHS = rbind(
+ XPy,
+ ZPy,
+ ZPy
+
+ )
>
> # 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] [,6]
[1,] 3.362958e-01 -2.302930e-01 6.971515e-02 2.675149e-02 1.060027e-01 9.646664e-02
[2,] -2.302930e-01 4.240019e-01 3.235525e-02 7.621335e-02 1.937088e-01 1.085686e-01
[3,] 6.971515e-02 3.235525e-02 3.338200e-01 -2.218598e-01 1.020704e-01 1.119602e-01
[4,] 2.675149e-02 7.621335e-02 -2.218598e-01 4.079596e-01 1.029648e-01 1.860998e-01
[5,] 1.060027e-01 1.937088e-01 1.020704e-01 1.029648e-01 2.997116e-01 2.050352e-01
[6,] 9.646664e-02 1.085686e-01 1.119602e-01 1.860998e-01 2.050352e-01 2.980600e-01
[7,] -3.221436e-02 -1.487315e-01 -9.972790e-02 -7.157314e-02 -1.809458e-01 -1.713010e-01
[8,] -1.212834e-01 -3.221494e-02 -8.531678e-02 -7.331967e-02 -1.534984e-01 -1.586365e-01
[9,] -8.459746e-02 -5.714882e-02 -5.305056e-02 -9.320243e-02 -1.417463e-01 -1.462530e-01
[10,] -1.549226e-01 -8.791262e-02 -1.683112e-01 -7.261134e-02 -2.428352e-01 -2.409226e-01
[11,] -1.460502e-01 -7.945741e-02 -7.871157e-02 -1.564158e-01 -2.255076e-01 -2.351274e-01
[12,] -7.019156e-02 -1.877674e-01 -1.450206e-01 -1.022485e-01 -2.579590e-01 -2.472691e-01
[13,] -1.612478e-01 -9.490402e-02 -1.542035e-01 -1.029546e-01 -2.561518e-01 -2.571581e-01
[14,] -6.971241e-02 -1.994368e-01 -9.510386e-02 -1.704180e-01 -2.691492e-01 -2.655219e-01
[15,] 1.658062e-17 -2.812342e-18 1.023140e-18 2.977664e-17 -2.177293e-17 -2.357044e-17
[16,] 6.812843e-18 4.417607e-17 -4.178285e-17 4.640101e-17 9.907072e-17 1.286827e-16
[17,] 6.271884e-17 2.579012e-17 5.588952e-17 2.573204e-17 6.072600e-18 1.048149e-17
[18,] -6.891580e-02 1.649668e-02 -6.483306e-02 1.183069e-02 -5.241912e-02 -5.300237e-02
[19,] -6.309608e-02 2.473172e-03 2.204704e-02 -9.483325e-02 -6.062291e-02 -7.278621e-02
[20,] 2.272884e-02 -8.840761e-02 -6.696106e-02 1.409512e-02 -6.567877e-02 -5.286593e-02
[21,] -6.599648e-02 1.972347e-02 -7.134430e-02 2.583526e-02 -4.627301e-02 -4.550904e-02
[22,] 3.242238e-02 -9.314286e-02 3.823423e-02 -9.978497e-02 -6.072048e-02 -6.155074e-02
[,7] [,8] [,9] [,10] [,11] [,12]
[1,] -3.221436e-02 -1.212834e-01 -8.459746e-02 -1.549226e-01 -1.460502e-01 -7.019156e-02
[2,] -1.487315e-01 -3.221494e-02 -5.714882e-02 -8.791262e-02 -7.945741e-02 -1.877674e-01
[3,] -9.972790e-02 -8.531678e-02 -5.305056e-02 -1.683112e-01 -7.871157e-02 -1.450206e-01
[4,] -7.157314e-02 -7.331967e-02 -9.320243e-02 -7.261134e-02 -1.564158e-01 -1.022485e-01
[5,] -1.809458e-01 -1.534984e-01 -1.417463e-01 -2.428352e-01 -2.255076e-01 -2.579590e-01
[6,] -1.713010e-01 -1.586365e-01 -1.462530e-01 -2.409226e-01 -2.351274e-01 -2.472691e-01
[7,] 6.392168e-01 2.363105e-02 5.143787e-02 3.083020e-01 8.428748e-02 3.446969e-01
[8,] 2.363105e-02 6.791768e-01 1.147783e-02 3.377197e-01 3.239027e-01 1.851742e-01
[9,] 5.143787e-02 1.147783e-02 6.513700e-01 6.826402e-02 3.060956e-01 1.844147e-01
[10,] 3.083020e-01 3.377197e-01 6.826402e-02 6.094300e-01 2.542925e-01 3.059343e-01
[11,] 8.428748e-02 3.239027e-01 3.060956e-01 2.542925e-01 5.887043e-01 3.289698e-01
[12,] 3.446969e-01 1.851742e-01 1.844147e-01 3.059343e-01 3.289698e-01 6.175962e-01
[13,] 1.845548e-01 2.075013e-01 3.222297e-01 3.558104e-01 3.125844e-01 2.854138e-01
[14,] 3.769941e-01 1.414741e-01 1.958175e-01 3.293019e-01 2.400521e-01 3.379905e-01
[15,] 2.395790e-17 4.883271e-18 -2.708100e-17 1.803468e-17 -1.615002e-17 1.904230e-18
[16,] -6.080268e-17 -1.220990e-16 -6.613368e-17 -1.055715e-16 -1.424444e-16 -1.283217e-16
[17,] -4.945417e-17 -3.400154e-17 -4.537775e-17 -6.495460e-17 -4.728673e-17 -5.206048e-17
[18,] -3.055726e-02 -3.616256e-02 6.671982e-02 -9.504654e-02 4.080270e-02 2.505500e-02
[19,] 6.633679e-02 -3.255401e-02 -3.378277e-02 4.677847e-02 -9.560950e-02 4.552758e-03
[20,] -2.046634e-02 1.368881e-02 6.777533e-03 2.956453e-02 -9.031318e-03 -9.165603e-02
[21,] 2.655683e-02 2.393825e-02 -5.049508e-02 2.200864e-02 1.389874e-02 3.452602e-02
[22,] -4.187001e-02 3.108952e-02 1.078050e-02 -3.305105e-03 4.993938e-02 2.752225e-02
[,13] [,14] [,15] [,16] [,17] [,18]
[1,] -1.612478e-01 -6.971241e-02 -1.641941e-15 -2.325870e-16 -1.574980e-17 -6.891580e-02
[2,] -9.490402e-02 -1.994368e-01 -1.646013e-15 -1.898534e-16 -4.427988e-17 1.649668e-02
[3,] -1.542035e-01 -9.510386e-02 4.117312e-16 1.324522e-16 1.315660e-16 -6.483306e-02
[4,] -1.029546e-01 -1.704180e-01 4.702551e-16 3.848965e-17 7.029055e-18 1.183069e-02
[5,] -2.561518e-01 -2.691492e-01 1.566591e-15 2.572675e-16 3.622674e-17 -5.241912e-02
[6,] -2.571581e-01 -2.655219e-01 -5.109661e-16 -5.578429e-17 -7.509013e-17 -5.300237e-02
[7,] 1.845548e-01 3.769941e-01 1.598326e-16 -1.212212e-16 -1.196552e-17 -3.055726e-02
[8,] 2.075013e-01 1.414741e-01 3.106246e-17 -5.859214e-17 -3.470884e-17 -3.616256e-02
[9,] 3.222297e-01 1.958175e-01 -1.256668e-18 9.417605e-17 6.656983e-17 6.671982e-02
[10,] 3.558104e-01 3.293019e-01 1.158805e-16 -9.870506e-17 -5.028129e-17 -9.504654e-02
[11,] 3.125844e-01 2.400521e-01 3.746834e-17 3.671045e-18 3.963854e-17 4.080270e-02
[12,] 2.854138e-01 3.379905e-01 1.332151e-16 -8.981793e-17 2.475122e-17 2.505500e-02
[13,] 6.135308e-01 3.869710e-01 6.518561e-17 9.330102e-18 -5.458000e-18 2.703240e-02
[14,] 3.869710e-01 6.594414e-01 1.182485e-16 -5.351458e-17 -1.239952e-18 9.434650e-03
[15,] -1.674072e-18 1.360012e-17 4.285714e-01 -1.249001e-16 -2.602085e-17 -1.611237e-18
[16,] -9.498857e-17 -9.607831e-17 6.938894e-17 4.285714e-01 4.727121e-17 8.768364e-18
[17,] -7.410980e-17 -6.773272e-17 2.949030e-17 -4.380177e-17 4.285714e-01 8.989137e-19
[18,] 2.703240e-02 9.434650e-03 -2.121626e-17 3.104273e-17 6.352098e-18 3.298300e-01
[19,] 3.515618e-02 6.800707e-02 5.976514e-18 2.529698e-18 -1.307603e-17 1.866951e-02
[20,] 4.421305e-02 3.536092e-02 -4.069808e-17 1.819295e-17 -4.029427e-17 2.391874e-02
[21,] -9.191546e-02 -1.718182e-02 1.247080e-17 -4.051109e-17 -2.003227e-17 4.271666e-02
[22,] -1.448617e-02 -9.562082e-02 -3.932732e-17 1.867247e-17 2.510164e-17 1.343650e-02
[,19] [,20] [,21] [,22]
[1,] -6.309608e-02 2.272884e-02 -6.599648e-02 3.242238e-02
[2,] 2.473172e-03 -8.840761e-02 1.972347e-02 -9.314286e-02
[3,] 2.204704e-02 -6.696106e-02 -7.134430e-02 3.823423e-02
[4,] -9.483325e-02 1.409512e-02 2.583526e-02 -9.978497e-02
[5,] -6.062291e-02 -6.567877e-02 -4.627301e-02 -6.072048e-02
[6,] -7.278621e-02 -5.286593e-02 -4.550904e-02 -6.155074e-02
[7,] 6.633679e-02 -2.046634e-02 2.655683e-02 -4.187001e-02
[8,] -3.255401e-02 1.368881e-02 2.393825e-02 3.108952e-02
[9,] -3.378277e-02 6.777533e-03 -5.049508e-02 1.078050e-02
[10,] 4.677847e-02 2.956453e-02 2.200864e-02 -3.305105e-03
[11,] -9.560950e-02 -9.031318e-03 1.389874e-02 4.993938e-02
[12,] 4.552758e-03 -9.165603e-02 3.452602e-02 2.752225e-02
[13,] 3.515618e-02 4.421305e-02 -9.191546e-02 -1.448617e-02
[14,] 6.800707e-02 3.536092e-02 -1.718182e-02 -9.562082e-02
[15,] 1.146659e-17 1.795250e-17 1.323016e-17 -1.324914e-18
[16,] 6.516845e-18 1.538432e-17 -9.098847e-19 -4.103621e-17
[17,] -1.359386e-18 3.590251e-18 7.196728e-18 2.516784e-17
[18,] 1.866951e-02 2.391874e-02 4.271666e-02 1.343650e-02
[19,] 3.421286e-01 2.302693e-02 2.403364e-02 2.071270e-02
[20,] 2.302693e-02 3.362828e-01 1.715788e-02 2.818508e-02
[21,] 2.403364e-02 1.715788e-02 3.260662e-01 1.859699e-02
[22,] 2.071270e-02 2.818508e-02 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
>
> # solutions for fixed effects and animal effects
> # hys
> sol_hys = sol[1 : no_lvl_f1]
>
> # parity
> sol_parity = sol[(no_lvl_f1 + 1) : (no_lvl_f1 + no_lvl_f2)]
>
> # animal
> sol_animal = sol[(no_lvl_f1 + no_lvl_f2 + 1) : (no_lvl_f1 + no_lvl_f2 + no_animals)]
>
> # permanent environmental
> sol_pe = sol[(no_lvl_f1 + no_lvl_f2 + no_animals + 1) : (no_lvl_f1 + no_lvl_f2 + no_animals * 2)]
>
> # print
> sol_hys
[1] 43.80219 87.86764 80.62658 80.63976
> sol_parity
[1] 131.6698 161.2663
> sol_animal
[1] 10.147570 -3.084153 -7.063418 13.580743 -18.206972 -18.386782 9.328432 24.193618
> sol_pe
[1] -1.468082e-15 1.496534e-14 -1.873529e-15 8.416979e+00 -7.145576e+00 -1.722850e+01
[7] -1.389647e+00 1.734674e+01
>
> # reliability(r2), accuracy(r), standard error of prediction(SEP)
>
> # diagonal elements of the generalized inverse of LHS for animal equation
> d_ani = diag(gi_LHS[(no_lvl_f1 + no_lvl_f2 + 1) : (no_lvl_f1 + no_lvl_f2 + no_animals),
+ (no_lvl_f1 + no_lvl_f2 + 1) : (no_lvl_f1 + no_lvl_f2 + no_animals)])
>
> # print
> d_ani
[1] 0.6392168 0.6791768 0.6513700 0.6094300 0.5887043 0.6175962 0.6135308 0.6594414
>
> # reliability
> rel = 1 - d_ani * alpha
>
> # print
> rel
[1] 0.10509649 0.04915244 0.08808198 0.14679806 0.17581391 0.13536537 0.14105693 0.07678208
>
> # accuracy
> acc = sqrt(rel)
>
> # print
> acc
[1] 0.3241859 0.2217035 0.2967861 0.3831423 0.4193017 0.3679203 0.3755755 0.2770958
>
> # standard error of prediction(SEP)
> sep = sqrt(d_ani * sigma_e)
>
> # 확인
> sep
[1] 4.230611 4.360843 4.270639 4.130864 4.060015 4.158448 4.144739 4.297017
>
> # partitioning of breeding values
>
> # yield deviation
> yd1 = ginv(ZPZ) %*% t(Z) %*% (y - X1 %*% sol_hys - X2 %*% sol_parity - Z %*% sol_pe)
>
> # print
> 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(ZPZ)
>
> # print
> a2
[1] 0 0 0 2 2 2 2 2
>
>
> # Parents average, progeny contribution
>
> # parents avearge
> pa1 = rep(0, no_animals)
>
> # progeny contribution numerator
> pc0 = rep(0, no_animals)
>
> # numerator of n3, denominator of progeny contribution
> a3 = rep(0, no_animals)
>
> # numerator of n1
> a1 = rep(0, 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 * alpha
+
+ } else if (sire != 0 & dam == 0) {
+ # 개체의 부만 알고 있을 경우
+
+ # PA
+ a1[i] = 4/3 * alpha
+ pa1[i] = sol_animal[sire]/2
+
+ # PC for sire
+ a3[sire] = a3[sire] + 0.5 * alpha * (2/3)
+ pc0[sire] = pc0[sire] + 0.5 * alpha * (2/3) * (2 * sol_animal[i])
+
+ } else if (sire == 0 & dam != 0) {
+ # dam known
+
+ # PA
+ a1[i] = 4/3 * alpha
+ pa1[i] = sol_animal[dam]/2
+
+ # PC for dam
+ a3[dam] = a3[dam] + 0.5 * alpha * (2/3)
+ pc0[dam] = pc0[dam] + 0.5 * alpha * (2/3) * (2 * sol_animal[i])
+
+ } else {
+ # both parents known
+
+ # PA
+ a1[i] = 2 * alpha
+ pa1[i] = (sol_animal[sire] + sol_animal[dam])/2
+
+ # PC for sire
+ a3[sire] = a3[sire] + 0.5 * alpha
+ pc0[sire] = pc0[sire] + 0.5 * alpha * (2 * sol_animal[i] - sol_animal[dam])
+
+ # PC for dam
+ a3[dam] = a3[dam] + 0.5 * alpha
+ pc0[dam] = pc0[dam] + 0.5 * alpha * (2 * sol_animal[i] - sol_animal[sire])
+
+ }
+ }
>
> # print
> 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 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 0.000000
>
> # denominator of n1, n2, n3, diagonal of animals in LHS
> n_de = a1 + a2 + a3
>
> # print
> 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)
>
> # print
> pa
[1] 0.000000 0.000000 0.000000 1.797961 -2.583018 -2.350659 1.658956 5.680501
>
> # yield deviation fraction of breeding values
> yd = yd1 * (a2 / n_de)
>
> # print
> 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 0.000000
>
> # progeny contribution fraction of breeding values
> pc = pc1 * (a3 / n_de)
>
> # print
> pc
[1] 10.147570 -3.084153 -7.063418 3.273490 -5.971781 0.000000 4.866867 0.000000
>
> # Progeny(Daughter) Yield Deviation(PYD, DYD)
>
> # n2 of progeny
> n2prog = a2 / (a1 + a2)
>
> # print
> n2prog
[1] 0.0000000 0.0000000 0.0000000 0.4166667 0.4166667 0.4166667 0.4166667 0.4166667
>
> # numerator of dyd : summation of u of progeny * n2 of progeny * (2 * YD - bv of mate)
> dyd_n = rep(0, no_animals)
> # denominator of dyd : summation of u of progeny * n2 of progeny
> dyd_d = rep(0, 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
+
+ } 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]
+
+ }
+ }
>
> # print
> dyd_n
[1] 29.438905 -3.904169 -20.070496 9.365739 -36.300400 0.000000 32.798081 0.000000
> dyd_d
[1] 1.2500000 0.8333333 0.8333333 0.4166667 0.4166667 0.0000000 0.4166667 0.0000000
>
> # dyd
> dyd = dyd_n / dyd_d
> dyd[is.nan(dyd) == TRUE] = 0
>
> # print
> dyd
[1] 23.551124 -4.685002 -24.084595 22.477774 -87.120960 0.000000 78.715394 0.000000
>
> # breeding values and fractions
> repeat_result = data.frame(animal = pedi[,1], animal_bv = sol_animal, rel = rel, acc = acc, sep = sep, pa = pa, yd = yd, pc = pc, sum_of_fr = pa + yd + pc, dyd = dyd)
>
> # print
> repeat_result
animal animal_bv rel acc sep pa yd pc sum_of_fr
1 1 10.147570 0.10509649 0.3241859 4.230611 0.000000 0.000000 10.147570 10.147570
2 2 -3.084153 0.04915244 0.2217035 4.360843 0.000000 0.000000 -3.084153 -3.084153
3 3 -7.063418 0.08808198 0.2967861 4.270639 0.000000 0.000000 -7.063418 -7.063418
4 4 13.580743 0.14679806 0.3831423 4.130864 1.797961 8.509292 3.273490 13.580743
5 5 -18.206972 0.17581391 0.4193017 4.060015 -2.583018 -9.652174 -5.971781 -18.206972
6 6 -18.386782 0.13536537 0.3679203 4.158448 -2.350659 -16.036123 0.000000 -18.386782
7 7 9.328432 0.14105693 0.3755755 4.144739 1.658956 2.802610 4.866867 9.328432
8 8 24.193618 0.07678208 0.2770958 4.297017 5.680501 18.513118 0.000000 24.193618
dyd
1 23.551124
2 -4.685002
3 -24.084595
4 22.477774
5 -87.120960
6 0.000000
7 78.715394
8 0.000000
>
> # 파일로 출력, 분리자는 ",", 따옴표 없고
> output_filename = gsub("[ -]", "", paste("repeat_result_", Sys.Date(), ".csv"))
> write.table(repeat_result, output_filename, sep=",", row.names = FALSE, quote = FALSE)
>
다음은 출력 파일을 엑셀에서 표시한 결과이다.
