돼지에서 어미의 한배새끼들이 어미로부터 같은 육성환경을 제공받는다면 공통 환경(common environment)을 제공받는다고 할 수 있다. 이러한 새끼들은 유전적인 원인 이외에 이러한 공통 환경 때문에 유사성이 발생하고 정확한 개체의 유전능력을 추정하기 위해서는 공통 환경 효과를 모형에 넣어야 한다. 어미의 육성 환경 중 유전적인 부분은 모체 효과(maternal effect)라고 하며, 공통 환경 효과의 특별 케이스라고 할 수 있다. 모체 효과에 대해서는 후에 다시 다루게 되고 여기서는 어미의 비유전적 환경이라고 할 수 있는 공통 환경 효과가 있는 모형의 육종가 구하기 프로그램이다.
자료는 다음과 같다.
6 2 1 90
7 2 2 70
8 2 2 65
9 4 2 98
10 4 1 106
11 4 2 60
12 4 2 80
13 5 1 100
14 5 2 85
15 5 1 68
개체, 어미, 성, 이유시체중
ce_data.txt로 저장
혈통은 다음과 같다.
1 0 0
2 0 0
3 0 0
4 0 0
5 0 0
6 1 2
7 1 2
8 1 2
9 3 4
10 3 4
11 3 4
12 3 4
13 1 5
14 1 5
15 1 5
개체, 아비, 어미
ce_pedi.txt 로 저장
육종가를 구하고, 신뢰도, 정확도, SEP를 구하고, 육종가를 분할하고, DYD를 구하는 R 프로그램은 다음과 같다.
# Common Environmental Effect Model - Date : 2020-07-01
# Linear Models for the Prediction of Animal Breeding Values, 3rd Edtion
# Raphael Mrode
# Example 4.2 p67
# MASS packages
if (!("MASS" %in% installed.packages())){
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/03_common_environment_R")
# print working directory
getwd()
# prdigree and data
# read pedigree : animal sire dam
pedi = read.table("ce_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("ce_data.txt", header = FALSE, sep = "", col.names = c("animal", "dam", "sex", "ww"), stringsAsFactors = FALSE)
# print
data
# variance component and ratio
sigma_a = 20
sigma_e = 65
sigma_c = 15
alpha = sigma_e/sigma_a
alpha2 = sigma_e/sigma_c
# print
sigma_a
sigma_e
sigma_c
alpha
alpha2
# design matrix
# design matrix of sex fixed effect
X = desgn(data[, 3])
# print
X
# number of levels of sex fixed levels
no_lvl_f1 = ncol(X)
# print
no_lvl_f1
# desing matrix of animal effect
Z = desgn(data[, 1])
# print
Z
# number of animals
no_animals = ncol(Z)
# print
no_animals
# desing matrix of common envrionmental effect(dam)
W = desgn(data[, 2])[,unique(data[,2])]
# print
W
# number of animals
no_ce = ncol(W)
# print
no_ce
# 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
XPW = t(X) %*% W
ZPX = t(Z) %*% X
ZPZ = t(Z) %*% Z
ZPW = t(Z) %*% W
WPX = t(W) %*% X
WPZ = t(W) %*% Z
WPW = t(W) %*% W
#print
XPX
XPZ
XPW
ZPX
ZPZ
ZPW
WPX
WPZ
WPW
LHS = rbind(
cbind(XPX, XPZ, XPW),
cbind(ZPX, ZPZ + ai * alpha, ZPW),
cbind(WPX, WPZ, WPW + diag(no_ce) * alpha2)
)
# print
LHS
# RHS construction
XPy = t(X) %*% y
ZPy = t(Z) %*% y
WPy = t(W) %*% y
# print
XPy
ZPy
WPy
RHS = rbind(
XPy,
ZPy,
WPy
)
# 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
# sex
sol_sex = sol[1 : no_lvl_f1]
# animal
sol_animal = sol[(no_lvl_f1 + 1) : (no_lvl_f1 + no_animals)]
# common environmental
sol_ce = sol[(no_lvl_f1 + no_animals + 1) : (no_lvl_f1 + no_animals + no_ce)]
# print
sol_sex
sol_animal
sol_ce
# 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 + 1) : (no_lvl_f1 + no_animals),
(no_lvl_f1 + 1) : (no_lvl_f1 + 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 - X %*% sol_sex - W %*% sol_ce)
# 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
common_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
common_result
# 파일로 출력, 분리자는 ",", 따옴표 없고
output_filename = gsub("[ -]", "", paste("common_result_", Sys.Date(), ".csv"))
write.table(common_result, output_filename, sep=",", row.names = FALSE, quote = FALSE)
다음은 실행 결과다.
> # Common Environmental Effect Model - Date : 2020-07-01
>
> # Linear Models for the Prediction of Animal Breeding Values, 3rd Edtion
>
> # Raphael Mrode
>
> # Example 4.2 p67
>
> # MASS packages
> if (!("MASS" %in% installed.packages())){
+ 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/03_common_environment_R")
>
> # print working directory
> getwd()
[1] "D:/temp/03_common_environment_R"
>
> # prdigree and data
>
> # read pedigree : animal sire dam
> pedi = read.table("ce_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 0 0
5 5 0 0
6 6 1 2
7 7 1 2
8 8 1 2
9 9 3 4
10 10 3 4
11 11 3 4
12 12 3 4
13 13 1 5
14 14 1 5
15 15 1 5
>
> # read data : animal, dam, sex, weaning_weight
> data = read.table("ce_data.txt", header = FALSE, sep = "", col.names = c("animal", "dam", "sex", "ww"), stringsAsFactors = FALSE)
>
> # print
> data
animal dam sex ww
1 6 2 1 90
2 7 2 2 70
3 8 2 2 65
4 9 4 2 98
5 10 4 1 106
6 11 4 2 60
7 12 4 2 80
8 13 5 1 100
9 14 5 2 85
10 15 5 1 68
>
> # variance component and ratio
> sigma_a = 20
> sigma_e = 65
> sigma_c = 15
>
> alpha = sigma_e/sigma_a
> alpha2 = sigma_e/sigma_c
>
> # print
> sigma_a
[1] 20
> sigma_e
[1] 65
> sigma_c
[1] 15
> alpha
[1] 3.25
> alpha2
[1] 4.333333
>
>
> # design matrix
>
> # design matrix of sex fixed effect
> X = desgn(data[, 3])
>
> # print
> X
[,1] [,2]
[1,] 1 0
[2,] 0 1
[3,] 0 1
[4,] 0 1
[5,] 1 0
[6,] 0 1
[7,] 0 1
[8,] 1 0
[9,] 0 1
[10,] 1 0
>
> # number of levels of sex fixed levels
> no_lvl_f1 = ncol(X)
>
> # print
> no_lvl_f1
[1] 2
>
> # desing matrix of animal effect
> Z = desgn(data[, 1])
>
> # print
> Z
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
[1,] 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
[2,] 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
[3,] 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
[4,] 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
[5,] 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
[6,] 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 1 0 0 0
[8,] 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 1 0
[10,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
>
> # number of animals
> no_animals = ncol(Z)
>
> # print
> no_animals
[1] 15
>
> # desing matrix of common envrionmental effect(dam)
> W = desgn(data[, 2])[,unique(data[,2])]
>
> # print
> W
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 1 0 0
[3,] 1 0 0
[4,] 0 1 0
[5,] 0 1 0
[6,] 0 1 0
[7,] 0 1 0
[8,] 0 0 1
[9,] 0 0 1
[10,] 0 0 1
>
> # number of animals
> no_ce = ncol(W)
>
> # print
> no_ce
[1] 3
>
> # observation
> y = data[, 4]
>
> # print
> y
[1] 90 70 65 98 106 60 80 100 85 68
>
> # inverse matrix of NRM
> ai = ainv(pedi)
>
> # print
> ai
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
[1,] 4.0 1.5 0 0 1.5 -1 -1 -1 0 0 0 0 -1 -1 -1
[2,] 1.5 2.5 0 0 0.0 -1 -1 -1 0 0 0 0 0 0 0
[3,] 0.0 0.0 3 2 0.0 0 0 0 -1 -1 -1 -1 0 0 0
[4,] 0.0 0.0 2 3 0.0 0 0 0 -1 -1 -1 -1 0 0 0
[5,] 1.5 0.0 0 0 2.5 0 0 0 0 0 0 0 -1 -1 -1
[6,] -1.0 -1.0 0 0 0.0 2 0 0 0 0 0 0 0 0 0
[7,] -1.0 -1.0 0 0 0.0 0 2 0 0 0 0 0 0 0 0
[8,] -1.0 -1.0 0 0 0.0 0 0 2 0 0 0 0 0 0 0
[9,] 0.0 0.0 -1 -1 0.0 0 0 0 2 0 0 0 0 0 0
[10,] 0.0 0.0 -1 -1 0.0 0 0 0 0 2 0 0 0 0 0
[11,] 0.0 0.0 -1 -1 0.0 0 0 0 0 0 2 0 0 0 0
[12,] 0.0 0.0 -1 -1 0.0 0 0 0 0 0 0 2 0 0 0
[13,] -1.0 0.0 0 0 -1.0 0 0 0 0 0 0 0 2 0 0
[14,] -1.0 0.0 0 0 -1.0 0 0 0 0 0 0 0 0 2 0
[15,] -1.0 0.0 0 0 -1.0 0 0 0 0 0 0 0 0 0 2
>
> # LHS, RHS
>
> # LHS construction
> XPX = t(X) %*% X
> XPZ = t(X) %*% Z
> XPW = t(X) %*% W
> ZPX = t(Z) %*% X
> ZPZ = t(Z) %*% Z
> ZPW = t(Z) %*% W
> WPX = t(W) %*% X
> WPZ = t(W) %*% Z
> WPW = t(W) %*% W
>
> #print
> XPX
[,1] [,2]
[1,] 4 0
[2,] 0 6
> XPZ
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
[1,] 0 0 0 0 0 1 0 0 0 1 0 0 1 0 1
[2,] 0 0 0 0 0 0 1 1 1 0 1 1 0 1 0
> XPW
[,1] [,2] [,3]
[1,] 1 1 2
[2,] 2 3 1
> ZPX
[,1] [,2]
[1,] 0 0
[2,] 0 0
[3,] 0 0
[4,] 0 0
[5,] 0 0
[6,] 1 0
[7,] 0 1
[8,] 0 1
[9,] 0 1
[10,] 1 0
[11,] 0 1
[12,] 0 1
[13,] 1 0
[14,] 0 1
[15,] 1 0
> ZPZ
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
[1,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[2,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[3,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[4,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[5,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[6,] 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
[7,] 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
[8,] 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
[9,] 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
[10,] 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
[11,] 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
[12,] 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
[13,] 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
[14,] 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
[15,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
> ZPW
[,1] [,2] [,3]
[1,] 0 0 0
[2,] 0 0 0
[3,] 0 0 0
[4,] 0 0 0
[5,] 0 0 0
[6,] 1 0 0
[7,] 1 0 0
[8,] 1 0 0
[9,] 0 1 0
[10,] 0 1 0
[11,] 0 1 0
[12,] 0 1 0
[13,] 0 0 1
[14,] 0 0 1
[15,] 0 0 1
> WPX
[,1] [,2]
[1,] 1 2
[2,] 1 3
[3,] 2 1
> WPZ
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
[1,] 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0
[2,] 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0
[3,] 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1
> WPW
[,1] [,2] [,3]
[1,] 3 0 0
[2,] 0 4 0
[3,] 0 0 3
>
> LHS = rbind(
+ cbind(XPX, XPZ, XPW),
+ cbind(ZPX, ZPZ + ai * alpha, ZPW),
+ cbind(WPX, WPZ, WPW + diag(no_ce) * alpha2)
+ )
>
> # print
> LHS
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
[1,] 4 0 0.000 0.000 0.00 0.00 0.000 1.00 0.00 0.00 0.00 1.00 0.00 0.00 1.00
[2,] 0 6 0.000 0.000 0.00 0.00 0.000 0.00 1.00 1.00 1.00 0.00 1.00 1.00 0.00
[3,] 0 0 13.000 4.875 0.00 0.00 4.875 -3.25 -3.25 -3.25 0.00 0.00 0.00 0.00 -3.25
[4,] 0 0 4.875 8.125 0.00 0.00 0.000 -3.25 -3.25 -3.25 0.00 0.00 0.00 0.00 0.00
[5,] 0 0 0.000 0.000 9.75 6.50 0.000 0.00 0.00 0.00 -3.25 -3.25 -3.25 -3.25 0.00
[6,] 0 0 0.000 0.000 6.50 9.75 0.000 0.00 0.00 0.00 -3.25 -3.25 -3.25 -3.25 0.00
[7,] 0 0 4.875 0.000 0.00 0.00 8.125 0.00 0.00 0.00 0.00 0.00 0.00 0.00 -3.25
[8,] 1 0 -3.250 -3.250 0.00 0.00 0.000 7.50 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[9,] 0 1 -3.250 -3.250 0.00 0.00 0.000 0.00 7.50 0.00 0.00 0.00 0.00 0.00 0.00
[10,] 0 1 -3.250 -3.250 0.00 0.00 0.000 0.00 0.00 7.50 0.00 0.00 0.00 0.00 0.00
[11,] 0 1 0.000 0.000 -3.25 -3.25 0.000 0.00 0.00 0.00 7.50 0.00 0.00 0.00 0.00
[12,] 1 0 0.000 0.000 -3.25 -3.25 0.000 0.00 0.00 0.00 0.00 7.50 0.00 0.00 0.00
[13,] 0 1 0.000 0.000 -3.25 -3.25 0.000 0.00 0.00 0.00 0.00 0.00 7.50 0.00 0.00
[14,] 0 1 0.000 0.000 -3.25 -3.25 0.000 0.00 0.00 0.00 0.00 0.00 0.00 7.50 0.00
[15,] 1 0 -3.250 0.000 0.00 0.00 -3.250 0.00 0.00 0.00 0.00 0.00 0.00 0.00 7.50
[16,] 0 1 -3.250 0.000 0.00 0.00 -3.250 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[17,] 1 0 -3.250 0.000 0.00 0.00 -3.250 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
[18,] 1 2 0.000 0.000 0.00 0.00 0.000 1.00 1.00 1.00 0.00 0.00 0.00 0.00 0.00
[19,] 1 3 0.000 0.000 0.00 0.00 0.000 0.00 0.00 0.00 1.00 1.00 1.00 1.00 0.00
[20,] 2 1 0.000 0.000 0.00 0.00 0.000 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00
[,16] [,17] [,18] [,19] [,20]
[1,] 0.00 1.00 1.000000 1.000000 2.000000
[2,] 1.00 0.00 2.000000 3.000000 1.000000
[3,] -3.25 -3.25 0.000000 0.000000 0.000000
[4,] 0.00 0.00 0.000000 0.000000 0.000000
[5,] 0.00 0.00 0.000000 0.000000 0.000000
[6,] 0.00 0.00 0.000000 0.000000 0.000000
[7,] -3.25 -3.25 0.000000 0.000000 0.000000
[8,] 0.00 0.00 1.000000 0.000000 0.000000
[9,] 0.00 0.00 1.000000 0.000000 0.000000
[10,] 0.00 0.00 1.000000 0.000000 0.000000
[11,] 0.00 0.00 0.000000 1.000000 0.000000
[12,] 0.00 0.00 0.000000 1.000000 0.000000
[13,] 0.00 0.00 0.000000 1.000000 0.000000
[14,] 0.00 0.00 0.000000 1.000000 0.000000
[15,] 0.00 0.00 0.000000 0.000000 1.000000
[16,] 7.50 0.00 0.000000 0.000000 1.000000
[17,] 0.00 7.50 0.000000 0.000000 1.000000
[18,] 0.00 0.00 7.333333 0.000000 0.000000
[19,] 0.00 0.00 0.000000 8.333333 0.000000
[20,] 1.00 1.00 0.000000 0.000000 7.333333
>
> # RHS construction
> XPy = t(X) %*% y
> ZPy = t(Z) %*% y
> WPy = t(W) %*% y
>
> # print
> XPy
[,1]
[1,] 364
[2,] 458
> ZPy
[,1]
[1,] 0
[2,] 0
[3,] 0
[4,] 0
[5,] 0
[6,] 90
[7,] 70
[8,] 65
[9,] 98
[10,] 106
[11,] 60
[12,] 80
[13,] 100
[14,] 85
[15,] 68
> WPy
[,1]
[1,] 225
[2,] 344
[3,] 253
>
> RHS = rbind(
+ XPy,
+ ZPy,
+ WPy
+
+ )
>
> # print
> RHS
[,1]
[1,] 364
[2,] 458
[3,] 0
[4,] 0
[5,] 0
[6,] 0
[7,] 0
[8,] 90
[9,] 70
[10,] 65
[11,] 98
[12,] 106
[13,] 60
[14,] 80
[15,] 100
[16,] 85
[17,] 68
[18,] 225
[19,] 344
[20,] 253
>
> # Solutions
>
> # generalized inverse of LHS
> gi_LHS = ginv(LHS)
>
> # print
> gi_LHS
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 0.44292085 0.13589744 -0.105239688 -0.041248606 -0.048606466 -0.048606466 -0.063991081
[2,] 0.13589744 0.34358974 -0.087179487 -0.051282051 -0.066666667 -0.066666667 -0.035897436
[3,] -0.10523969 -0.08717949 0.286733556 -0.011148272 0.020958751 0.020958751 -0.009810479
[4,] -0.04124861 -0.05128205 -0.011148272 0.285395764 0.011148272 0.011148272 0.011148272
[5,] -0.04860647 -0.06666667 0.020958751 0.011148272 0.286733556 -0.020958751 0.009810479
[6,] -0.04860647 -0.06666667 0.020958751 0.011148272 -0.020958751 0.286733556 0.009810479
[7,] -0.06399108 -0.03589744 -0.009810479 0.011148272 0.009810479 0.009810479 0.286733556
[8,] -0.11428465 -0.06786325 0.135681903 0.128799703 0.018164251 0.018164251 0.006882200
[9,] -0.07334820 -0.09555556 0.133273876 0.130137495 0.020572278 0.020572278 0.003136381
[10,] -0.07334820 -0.09555556 0.133273876 0.130137495 0.020572278 0.020572278 0.003136381
[11,] -0.05052397 -0.09025641 0.025596433 0.014269788 0.128249721 0.128249721 0.011326644
[12,] -0.09146042 -0.06256410 0.028004459 0.012931996 0.125841695 0.125841695 0.015072464
[13,] -0.05052397 -0.09025641 0.025596433 0.014269788 0.128249721 0.128249721 0.011326644
[14,] -0.05052397 -0.09025641 0.025596433 0.014269788 0.128249721 0.128249721 0.011326644
[15,] -0.11959123 -0.06427350 0.135994054 0.003270160 0.017852100 0.017852100 0.132723894
[16,] -0.07865478 -0.09196581 0.133586027 0.004607952 0.020260126 0.020260126 0.128978075
[17,] -0.11959123 -0.06427350 0.135994054 0.003270160 0.017852100 0.017852100 0.132723894
[18,] -0.06187291 -0.07692308 -0.016722408 -0.033444816 0.016722408 0.016722408 0.016722408
[19,] -0.07290970 -0.10000000 0.031438127 0.016722408 -0.031438127 -0.031438127 0.014715719
[20,] -0.09598662 -0.05384615 -0.014715719 0.016722408 0.014715719 0.014715719 -0.031438127
[,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] -0.11428465 -0.073348198 -0.073348198 -0.05052397 -0.09146042 -0.05052397 -0.05052397
[2,] -0.06786325 -0.095555556 -0.095555556 -0.09025641 -0.06256410 -0.09025641 -0.09025641
[3,] 0.13568190 0.133273876 0.133273876 0.02559643 0.02800446 0.02559643 0.02559643
[4,] 0.12879970 0.130137495 0.130137495 0.01426979 0.01293200 0.01426979 0.01426979
[5,] 0.01816425 0.020572278 0.020572278 0.12824972 0.12584169 0.12824972 0.12824972
[6,] 0.01816425 0.020572278 0.020572278 0.12824972 0.12584169 0.12824972 0.12824972
[7,] 0.00688220 0.003136381 0.003136381 0.01132664 0.01507246 0.01132664 0.01132664
[8,] 0.26818927 0.128666419 0.128666419 0.02115793 0.02734745 0.02115793 0.02115793
[9,] 0.12866642 0.264960733 0.131627400 0.02645559 0.02349461 0.02645559 0.02645559
[10,] 0.12866642 0.131627400 0.264960733 0.02645559 0.02349461 0.02645559 0.02645559
[11,] 0.02115793 0.026455593 0.026455593 0.26163657 0.12300557 0.12830323 0.12830323
[12,] 0.02734745 0.023494612 0.023494612 0.12300557 0.26019175 0.12300557 0.12300557
[13,] 0.02115793 0.026455593 0.026455593 0.12830323 0.12300557 0.26163657 0.12830323
[14,] 0.02115793 0.026455593 0.026455593 0.12830323 0.12300557 0.12830323 0.26163657
[15,] 0.07563929 0.068263595 0.068263595 0.02047120 0.02784690 0.02047120 0.02047120
[16,] 0.06944977 0.071224576 0.071224576 0.02576886 0.02399405 0.02576886 0.02576886
[17,] 0.07563929 0.068263595 0.068263595 0.02047120 0.02784690 0.02047120 0.02047120
[18,] -0.03756968 -0.035562988 -0.035562988 0.02140468 0.01939799 0.02140468 0.02140468
[19,] 0.02724638 0.030858417 0.030858417 -0.03839465 -0.04200669 -0.03839465 -0.03839465
[20,] 0.01032330 0.004704571 0.004704571 0.01698997 0.02260870 0.01698997 0.01698997
[,15] [,16] [,17] [,18] [,19] [,20]
[1,] -0.11959123 -0.078654775 -0.11959123 -0.061872910 -0.07290970 -0.095986622
[2,] -0.06427350 -0.091965812 -0.06427350 -0.076923077 -0.10000000 -0.053846154
[3,] 0.13599405 0.133586027 0.13599405 -0.016722408 0.03143813 -0.014715719
[4,] 0.00327016 0.004607952 0.00327016 -0.033444816 0.01672241 0.016722408
[5,] 0.01785210 0.020260126 0.01785210 0.016722408 -0.03143813 0.014715719
[6,] 0.01785210 0.020260126 0.01785210 0.016722408 -0.03143813 0.014715719
[7,] 0.13272389 0.128978075 0.13272389 0.016722408 0.01471572 -0.031438127
[8,] 0.07563929 0.069449771 0.07563929 -0.037569677 0.02724638 0.010323300
[9,] 0.06826359 0.071224576 0.06826359 -0.035562988 0.03085842 0.004704571
[10,] 0.06826359 0.071224576 0.06826359 -0.035562988 0.03085842 0.004704571
[11,] 0.02047120 0.025768859 0.02047120 0.021404682 -0.03839465 0.016989967
[12,] 0.02784690 0.023994054 0.02784690 0.019397993 -0.04200669 0.022608696
[13,] 0.02047120 0.025768859 0.02047120 0.021404682 -0.03839465 0.016989967
[14,] 0.02047120 0.025768859 0.02047120 0.021404682 -0.03839465 0.016989967
[15,] 0.26994773 0.129238697 0.13661439 0.004905240 0.02677815 -0.031683389
[16,] 0.12923870 0.264346835 0.12923870 0.006911929 0.03039019 -0.037302118
[17,] 0.13661439 0.129238697 0.26994773 0.004905240 0.02677815 -0.031683389
[18,] 0.00490524 0.006911929 0.00490524 0.180602007 0.02508361 0.025083612
[19,] 0.02677815 0.030390190 0.02677815 0.025083612 0.18361204 0.022073579
[20,] -0.03168339 -0.037302118 -0.03168339 0.025083612 0.02207358 0.183612040
>
> # solution
> sol = gi_LHS %*% RHS
>
> # print
> sol
[,1]
[1,] 91.4931401
[2,] 75.7644444
[3,] -1.4407729
[4,] -1.1748792
[5,] 1.4407729
[6,] 1.4407729
[7,] -0.2658937
[8,] -1.0975588
[9,] -1.6670660
[10,] -2.3337327
[11,] 3.9252560
[12,] 2.8947633
[13,] -1.1414106
[14,] 1.5252560
[15,] 0.4478712
[16,] 0.5450306
[17,] -3.8187955
[18,] -1.7623188
[19,] 2.1611594
[20,] -0.3988406
>
> # solutions for fixed effects and animal effects
> # sex
> sol_sex = sol[1 : no_lvl_f1]
>
> # animal
> sol_animal = sol[(no_lvl_f1 + 1) : (no_lvl_f1 + no_animals)]
>
> # common environmental
> sol_ce = sol[(no_lvl_f1 + no_animals + 1) : (no_lvl_f1 + no_animals + no_ce)]
>
> # print
> sol_sex
[1] 91.49314 75.76444
> sol_animal
[1] -1.4407729 -1.1748792 1.4407729 1.4407729 -0.2658937 -1.0975588 -1.6670660 -2.3337327
[9] 3.9252560 2.8947633 -1.1414106 1.5252560 0.4478712 0.5450306 -3.8187955
> sol_ce
[1] -1.7623188 2.1611594 -0.3988406
>
> # 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 + 1) : (no_lvl_f1 + no_animals),
+ (no_lvl_f1 + 1) : (no_lvl_f1 + no_animals)])
>
> # print
> d_ani
[1] 0.2867336 0.2853958 0.2867336 0.2867336 0.2867336 0.2681893 0.2649607 0.2649607 0.2616366
[10] 0.2601918 0.2616366 0.2616366 0.2699477 0.2643468 0.2699477
>
> # reliability
> rel = 1 - d_ani * alpha
>
> # print
> rel
[1] 0.06811594 0.07246377 0.06811594 0.06811594 0.06811594 0.12838486 0.13887762 0.13887762
[9] 0.14968116 0.15437681 0.14968116 0.14968116 0.12266989 0.14087279 0.12266989
>
> # accuracy
> acc = sqrt(rel)
>
> # print
> acc
[1] 0.2609903 0.2691910 0.2609903 0.2609903 0.2609903 0.3583083 0.3726629 0.3726629 0.3868865
[10] 0.3929081 0.3868865 0.3868865 0.3502426 0.3753302 0.3502426
>
> # standard error of prediction(SEP)
> sep = sqrt(d_ani * sigma_e)
>
> # 확인
> sep
[1] 4.317138 4.307055 4.317138 4.317138 4.317138 4.175201 4.149994 4.149994 4.123879 4.112477
[11] 4.123879 4.123879 4.188866 4.145183 4.188866
>
> # partitioning of breeding values
>
> # yield deviation
> yd1 = ginv(ZPZ) %*% t(Z) %*% (y - X %*% sol_sex - W %*% sol_ce)
>
> # print
> yd1
[,1]
[1,] 0.0000000
[2,] 0.0000000
[3,] 0.0000000
[4,] 0.0000000
[5,] 0.0000000
[6,] 0.2691787
[7,] -4.0021256
[8,] -9.0021256
[9,] 20.0743961
[10,] 12.3457005
[11,] -17.9256039
[12,] 2.0743961
[13,] 8.9057005
[14,] 9.6343961
[15,] -23.0942995
>
> # numerator of n2
> a2 = diag(ZPZ)
>
> # print
> a2
[1] 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1
>
>
> # 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] 3.25 3.25 3.25 3.25 3.25 6.50 6.50 6.50 6.50 6.50 6.50 6.50 6.50 6.50 6.50
> pa1
[1] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 -1.3078261 -1.3078261 -1.3078261
[9] 1.4407729 1.4407729 1.4407729 1.4407729 -0.8533333 -0.8533333 -0.8533333
> a3
[1] 9.750 4.875 6.500 6.500 4.875 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
> pc0
[1] -18.730048 -9.545894 14.047536 14.047536 -2.160386 0.000000 0.000000 0.000000
[9] 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
>
> # denominator of n1, n2, n3, diagonal of animals in LHS
> n_de = a1 + a2 + a3
>
> # print
> n_de
[1] 13.000 8.125 9.750 9.750 8.125 7.500 7.500 7.500 7.500 7.500 7.500 7.500 7.500
[14] 7.500 7.500
>
> # parents average fraction of breeding values
> pa = pa1 * (a1 / n_de)
>
> # print
> pa
[1] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 -1.1334493 -1.1334493 -1.1334493
[9] 1.2486699 1.2486699 1.2486699 1.2486699 -0.7395556 -0.7395556 -0.7395556
>
> # yield deviation fraction of breeding values
> yd = yd1 * (a2 / n_de)
>
> # print
> yd
[,1]
[1,] 0.0000000
[2,] 0.0000000
[3,] 0.0000000
[4,] 0.0000000
[5,] 0.0000000
[6,] 0.0358905
[7,] -0.5336167
[8,] -1.2002834
[9,] 2.6765862
[10,] 1.6460934
[11,] -2.3900805
[12,] 0.2765862
[13,] 1.1874267
[14,] 1.2845862
[15,] -3.0792399
>
> # progeny contribution
> pc1 = pc0 / a3
> pc1[is.nan(pc1) == TRUE] = 0
> pc1
[1] -1.9210306 -1.9581320 2.1611594 2.1611594 -0.4431562 0.0000000 0.0000000 0.0000000
[9] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
>
> # progeny contribution fraction of breeding values
> pc = pc1 * (a3 / n_de)
>
> # print
> pc
[1] -1.4407729 -1.1748792 1.4407729 1.4407729 -0.2658937 0.0000000 0.0000000 0.0000000
[9] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
>
> # Progeny(Daughter) Yield Deviation(PYD, DYD)
>
> # n2 of progeny
> n2prog = a2 / (a1 + a2)
>
> # print
> n2prog
[1] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.1333333 0.1333333 0.1333333 0.1333333
[10] 0.1333333 0.1333333 0.1333333 0.1333333 0.1333333 0.1333333
>
> # 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] -4.0341643 -2.8197101 3.6499581 3.6499581 -0.6381449 0.0000000 0.0000000 0.0000000
[9] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
> dyd_d
[1] 0.8000000 0.4000000 0.5333333 0.5333333 0.4000000 0.0000000 0.0000000 0.0000000 0.0000000
[10] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
>
> # dyd
> dyd = dyd_n / dyd_d
> dyd[is.nan(dyd) == TRUE] = 0
>
> # print
> dyd
[1] -5.042705 -7.049275 6.843671 6.843671 -1.595362 0.000000 0.000000 0.000000 0.000000
[10] 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
>
> # breeding values and fractions
> common_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
> common_result
animal animal_bv rel acc sep pa yd pc sum_of_fr
1 1 -1.4407729 0.06811594 0.2609903 4.317138 0.0000000 0.0000000 -1.4407729 -1.4407729
2 2 -1.1748792 0.07246377 0.2691910 4.307055 0.0000000 0.0000000 -1.1748792 -1.1748792
3 3 1.4407729 0.06811594 0.2609903 4.317138 0.0000000 0.0000000 1.4407729 1.4407729
4 4 1.4407729 0.06811594 0.2609903 4.317138 0.0000000 0.0000000 1.4407729 1.4407729
5 5 -0.2658937 0.06811594 0.2609903 4.317138 0.0000000 0.0000000 -0.2658937 -0.2658937
6 6 -1.0975588 0.12838486 0.3583083 4.175201 -1.1334493 0.0358905 0.0000000 -1.0975588
7 7 -1.6670660 0.13887762 0.3726629 4.149994 -1.1334493 -0.5336167 0.0000000 -1.6670660
8 8 -2.3337327 0.13887762 0.3726629 4.149994 -1.1334493 -1.2002834 0.0000000 -2.3337327
9 9 3.9252560 0.14968116 0.3868865 4.123879 1.2486699 2.6765862 0.0000000 3.9252560
10 10 2.8947633 0.15437681 0.3929081 4.112477 1.2486699 1.6460934 0.0000000 2.8947633
11 11 -1.1414106 0.14968116 0.3868865 4.123879 1.2486699 -2.3900805 0.0000000 -1.1414106
12 12 1.5252560 0.14968116 0.3868865 4.123879 1.2486699 0.2765862 0.0000000 1.5252560
13 13 0.4478712 0.12266989 0.3502426 4.188866 -0.7395556 1.1874267 0.0000000 0.4478712
14 14 0.5450306 0.14087279 0.3753302 4.145183 -0.7395556 1.2845862 0.0000000 0.5450306
15 15 -3.8187955 0.12266989 0.3502426 4.188866 -0.7395556 -3.0792399 0.0000000 -3.8187955
dyd
1 -5.042705
2 -7.049275
3 6.843671
4 6.843671
5 -1.595362
6 0.000000
7 0.000000
8 0.000000
9 0.000000
10 0.000000
11 0.000000
12 0.000000
13 0.000000
14 0.000000
15 0.000000
>
> # 파일로 출력, 분리자는 ",", 따옴표 없고
> output_filename = gsub("[ -]", "", paste("common_result_", Sys.Date(), ".csv"))
> write.table(common_result, output_filename, sep=",", row.names = FALSE, quote = FALSE)
>
다음은 출력한 파일의 내용이다.
animal,animal_bv,rel,acc,sep,pa,yd,pc,sum_of_fr,dyd
1,-1.44077294685996,0.0681159420289843,0.260990310220484,4.3171380750933,0,0,-1.44077294685996,-1.44077294685996,-5.04270531400974
2,-1.17487922705307,0.0724637681159419,0.269190951029083,4.30705521646532,0,0,-1.17487922705312,-1.17487922705312,-7.04927536231885
3,1.44077294685983,0.068115942028986,0.260990310220487,4.3171380750933,0,0,1.44077294685986,1.44077294685986,6.84367149758441
4,1.44077294685994,0.0681159420289865,0.260990310220488,4.3171380750933,0,0,1.44077294685993,1.44077294685993,6.84367149758452
5,-0.2658937198067,0.0681159420289857,0.260990310220486,4.3171380750933,0,0,-0.26589371980672,-0.26589371980672,-1.59536231884048
6,-1.09755877616754,0.128384863123992,0.358308335270047,4.17520092181444,-1.13344927536231,0.0358904991948314,0,-1.09755877616748,0
7,-1.66706602254426,0.138877616747182,0.37266287277804,4.14999369458031,-1.13344927536231,-0.533616747181963,0,-1.66706602254428,0
8,-2.33373268921094,0.138877616747181,0.37266287277804,4.14999369458031,-1.13344927536231,-1.20028341384863,0,-2.33373268921094,0
9,3.92525603864735,0.14968115942029,0.386886494233503,4.12387885510646,1.24866988727857,2.67658615136876,0,3.92525603864732,0
10,2.89476328502408,0.154376811594203,0.392908146510356,4.11247659788064,1.24866988727857,1.64609339774555,0,2.89476328502412,0
11,-1.14141062801932,0.14968115942029,0.386886494233502,4.12387885510646,1.24866988727857,-2.39008051529791,0,-1.14141062801934,0
12,1.52525603864734,0.14968115942029,0.386886494233502,4.12387885510646,1.24866988727857,0.276586151368758,0,1.52525603864732,0
13,0.44787117552333,0.122669887278583,0.350242612025697,4.18886646414377,-0.739555555555552,1.1874267310789,0,0.447871175523351,0
14,0.545030595813275,0.140872785829307,0.375330235698255,4.14518326294675,-0.739555555555552,1.28458615136877,0,0.545030595813223,0
15,-3.81879549114334,0.122669887278583,0.350242612025697,4.18886646414377,-0.739555555555552,-3.07923993558776,0,-3.81879549114332,0
csv 파일이므로 엑셀에서 열면 다음과 같다.
'Animal Breeding > R for Genetic Evaluation' 카테고리의 다른 글
| R을 이용한 다형질 개체 모형(same model and no missing record)의 육종가 구하기 프로그래밍 (0) | 2020.07.03 |
|---|---|
| R을 이용하여 다형질 개체 모형(multiple trait animal model)에 대한 blupf90 분석 결과 정리하기 (0) | 2020.07.02 |
| R을 이용하여 공통 환경 효과 모형(common environmental effect model)에 대한 blupf90 분석 결과 정리하기 (0) | 2020.07.01 |
| R을 이용한 반복 모형(Repeatability Model)의 육종가 구하기 프로그래밍 (0) | 2020.06.14 |
| R을 이용하여 반복 모형(repeatability model)에 대한 blupf90 분석 결과 정리하기 (0) | 2020.06.13 |