R을 이용하여 단형질 개체 모형의 육종가, 신뢰도, 정확도, 육종가 분할 및 DYD(daughter yield deviation)을 구하는 프로그래밍을 한다. BLUPF90, ASREML, MIXBLUP, WOMBAT 등과 같은 전문적인 유전 능력 평가 프로그램을 이용할 경우 육종가(정확도 포함)를 구할 수 있으나, 그 과정을 이해하려면 R로 해야만 한다. 그리고 유전 능력 평가 프로그램도 육종가 분할이나 DYD를 구해 주진 않는다. 자료가 클 경우 본 프로그램을 이용하여 육종가를 구하는 것은 불가능한다. 그러나 육종가를 분할할 때 또는 DYD를 구해야 할 경우, 또는 알고리듬을 이해해야 하는 경우 유용할 것으로 생각한다.
-
혈통 자료 - stam_pedi.txt
animal, sire, dam
1, 0, 0
2, 0, 0
3, 0, 0
4, 1, 0
5, 3, 2
6, 1, 2
7, 4, 5
8, 3, 6
-
표현형 자료 - stam_data.txt
animal, sex, gain
4, 1, 4.5
5, 2, 2.9
6, 2, 3.9
7, 1, 3.5
8, 1, 5
-
프로그램
# Single Trait Animal Model - Date : 2020-06-01
# Linear Models for the Prediction of Animal Breeding Values, 3rd Edtion
# Raphael Mrode
# Example 3.1 p37
# 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/stam")
# print working directory
getwd()
# prdigree and data
# read pedigree : animal sire dam
pedi = read.table("stam_pedi.txt", header = TRUE, sep = ",", stringsAsFactors = FALSE)
# print
pedi
# read data : animal, sex, pre-weaning gain
data = read.table("stam_data.txt", header = TRUE, sep = ",", stringsAsFactors = FALSE)
# print
data
# variance component and ratio
sigma_a = 20
sigma_e = 40
alpha = sigma_e/sigma_a
# print
sigma_a
sigma_e
alpha
# design matrix
# design matrix of fixed effect
X = desgn(data[, 2])
# 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[, 3]
# print
y
# inverse matrix of NRM
ai = ainv(pedi)
# print
ai
# LHS, RHS
# LHS construction
LHS = rbind(
cbind(t(X) %*% X, t(X) %*% Z),
cbind(t(Z) %*% X, t(Z) %*% Z + ai * alpha))
# print
LHS
# RHS contruction
RHS = rbind(t(X) %*% y, t(Z) %*% y)
# 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_fixed = sol[1 : no_lvl_f]
sol_animal = sol[(no_lvl_f + 1) : (no_lvl_f + no_animals)]
# 확인
sol_fixed
sol_animal
# 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_f + 1) : (no_lvl_f + no_animals), (no_lvl_f + 1) : (no_lvl_f + 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(t(Z) %*% Z) %*% t(Z) %*% (y - X %*% sol_fixed)
# print
yd1
# numerator of n2
a2 = diag(t(Z) %*% Z)
# 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_nu = a1 + a2 + a3
# print
n_nu
# parents average fraction of breeding values
pa = pa1 * (a1 / n_nu)
# print
pa
# yield deviation fraction of breeding values
yd = yd1 * (a2 / n_nu)
# print
yd
# progeny contribution
pc1 = pc0 / a3
pc1[is.nan(pc1) == TRUE] = 0
pc1
# progeny contribution fraction of breeding values
pc = pc1 * (a3 / n_nu)
# 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
stam_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
stam_result
# 파일로 출력, 분리자는 ",", 따옴표 없고
output_filename = gsub("[ -]", "", paste("stam_result_", Sys.Date(), ".csv"))
write.table(stam_result, output_filename, sep=",", row.names = FALSE, quote = FALSE)
-
실행 결과
> # Single Trait Animal Model - Date : 2020-06-01
>
> # Linear Models for the Prediction of Animal Breeding Values, 3rd Edtion
>
> # Raphael Mrode
>
> # Example 3.1 p37
>
> # 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/stam")
>
> # print working directory
> getwd()
[1] "D:/temp/stam"
>
>
> # prdigree and data
>
> # read pedigree : animal sire dam
> pedi = read.table("stam_pedi.txt", header = TRUE, sep = ",", stringsAsFactors = FALSE)
>
> # 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, sex, pre-weaning gain
> data = read.table("stam_data.txt", header = TRUE, sep = ",", stringsAsFactors = FALSE)
>
> # print
> data
animal sex gain
1 4 1 4.5
2 5 2 2.9
3 6 2 3.9
4 7 1 3.5
5 8 1 5.0
>
> # variance component and ratio
> sigma_a = 20
> sigma_e = 40
> alpha = sigma_e/sigma_a
>
> # print
> sigma_a
[1] 20
> sigma_e
[1] 40
> alpha
[1] 2
>
> # design matrix
>
> # design matrix of fixed effect
> X = desgn(data[, 2])
>
> # print
> X
[,1] [,2]
[1,] 1 0
[2,] 0 1
[3,] 0 1
[4,] 1 0
[5,] 1 0
>
> # number of levels of fixed effect
> no_lvl_f = ncol(X)
>
> # print
> no_lvl_f
[1] 2
>
> # 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 0 1 0 0 0
[3,] 0 0 0 0 0 1 0 0
[4,] 0 0 0 0 0 0 1 0
[5,] 0 0 0 0 0 0 0 1
>
> # number of animals
> no_animals = ncol(Z)
>
> # print
> no_animals
[1] 8
>
>
> # observation
> y = data[, 3]
>
> # print
> y
[1] 4.5 2.9 3.9 3.5 5.0
>
> # 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
> LHS = rbind(
+ cbind(t(X) %*% X, t(X) %*% Z),
+ cbind(t(Z) %*% X, t(Z) %*% Z + ai * alpha))
>
> # print
> LHS
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 3 0 0.000000 0 0 1.000000 0 0 1 1
[2,] 0 2 0.000000 0 0 0.000000 1 1 0 0
[3,] 0 0 3.666667 1 0 -1.333333 0 -2 0 0
[4,] 0 0 1.000000 4 1 0.000000 -2 -2 0 0
[5,] 0 0 0.000000 1 4 0.000000 -2 1 0 -2
[6,] 1 0 -1.333333 0 0 4.666667 1 0 -2 0
[7,] 0 1 0.000000 -2 -2 1.000000 6 0 -2 0
[8,] 0 1 -2.000000 -2 1 0.000000 0 6 0 -2
[9,] 1 0 0.000000 0 0 -2.000000 -2 0 5 0
[10,] 1 0 0.000000 0 -2 0.000000 0 -2 0 5
>
> # RHS contruction
> RHS = rbind(t(X) %*% y, t(Z) %*% y)
>
> # print
> RHS
[,1]
[1,] 13.0
[2,] 6.8
[3,] 0.0
[4,] 0.0
[5,] 0.0
[6,] 4.5
[7,] 2.9
[8,] 3.9
[9,] 3.5
[10,] 5.0
>
> # Solutions
>
> # generalized inverse of LHS
> gi_LHS = ginv(LHS)
>
> # print
> gi_LHS
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 0.59556097 0.15730213 -0.164119413 -0.083624565 -0.13059083 -0.26455749 -0.14827804 -0.16632621 -0.2842464
[2,] 0.15730213 0.80245865 -0.132863260 -0.241250738 -0.11196840 -0.08730803 -0.29891465 -0.30600266 -0.1859495
[3,] -0.16411941 -0.13286326 0.471094211 0.006928037 0.03264668 0.21954371 0.04495225 0.22077427 0.1386223
[4,] -0.08362457 -0.24125074 0.006928037 0.492095721 -0.01030797 0.02039033 0.23734577 0.24515571 0.1198194
[5,] -0.13059083 -0.11196840 0.032646682 -0.010307967 0.45645878 0.04812709 0.20132326 0.02261354 0.1258983
[6,] -0.26455749 -0.08730803 0.219543709 0.020390333 0.04812709 0.42768015 0.04704420 0.12757186 0.2428012
[7,] -0.14827804 -0.29891465 0.044952254 0.237345770 0.20132326 0.04704420 0.42810675 0.16972255 0.2197160
[8,] -0.16632621 -0.30600266 0.220774267 0.245155707 0.02261354 0.12757186 0.16972255 0.44228277 0.1521830
[9,] -0.28424641 -0.18594950 0.138622268 0.119819354 0.12589831 0.24280124 0.21971599 0.15218301 0.4418562
[10,] -0.23787901 -0.19864885 0.134192262 0.110664009 0.21774710 0.12319108 0.17807393 0.21922376 0.1680818
[,10]
[1,] -0.2378790
[2,] -0.1986488
[3,] 0.1341923
[4,] 0.1106640
[5,] 0.2177471
[6,] 0.1231911
[7,] 0.1780739
[8,] 0.2192238
[9,] 0.1680818
[10,] 0.4223641
>
> # solution
> sol = gi_LHS %*% RHS
>
> # print
> sol
[,1]
[1,] 4.358502330
[2,] 3.404430006
[3,] 0.098444576
[4,] -0.018770099
[5,] -0.041084203
[6,] -0.008663123
[7,] -0.185732099
[8,] 0.176872088
[9,] -0.249458555
[10,] 0.182614688
>
> # solutions for fixed effects and animal effects
> sol_fixed = sol[1 : no_lvl_f]
> sol_animal = sol[(no_lvl_f + 1) : (no_lvl_f + no_animals)]
>
> # 확인
> sol_fixed
[1] 4.358502 3.404430
> sol_animal
[1] 0.098444576 -0.018770099 -0.041084203 -0.008663123 -0.185732099 0.176872088 -0.249458555 0.182614688
>
> # 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_f + 1) : (no_lvl_f + no_animals), (no_lvl_f + 1) : (no_lvl_f + no_animals)])
>
> # print
> d_ani
[1] 0.4710942 0.4920957 0.4564588 0.4276802 0.4281067 0.4422828 0.4418562 0.4223641
>
> # reliability
> rel = 1 - d_ani * alpha
>
> # print
> rel
[1] 0.05781158 0.01580856 0.08708243 0.14463969 0.14378651 0.11543447 0.11628766 0.15527171
>
> # accuracy
> acc = sqrt(rel)
>
> # print
> acc
[1] 0.2404404 0.1257321 0.2950973 0.3803153 0.3791919 0.3397565 0.3410098 0.3940453
>
> # standard error of prediction(SEP)
> sep = sqrt(d_ani * sigma_e)
>
> # 확인
> sep
[1] 4.340941 4.436646 4.272979 4.136086 4.138148 4.206104 4.204075 4.110300
>
> # partitioning of breeding values
>
> # yield deviation
> yd1 = ginv(t(Z) %*% Z) %*% t(Z) %*% (y - X %*% sol_fixed)
>
> # print
> yd1
[,1]
[1,] 0.0000000
[2,] 0.0000000
[3,] 0.0000000
[4,] 0.1414977
[5,] -0.5044300
[6,] 0.4955700
[7,] -0.8585023
[8,] 0.6414977
>
> # numerator of n2
> a2 = diag(t(Z) %*% Z)
>
> # print
> a2
[1] 0 0 0 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] 2.000000 2.000000 2.000000 2.666667 4.000000 4.000000 4.000000 4.000000
> pa1
[1] 0.00000000 0.00000000 0.00000000 0.04922229 -0.02992715 0.03983724 -0.09719761 0.06789394
> a3
[1] 1.666667 2.000000 2.000000 1.000000 1.000000 1.000000 0.000000 0.000000
> pc0
[1] 0.3609634 -0.0750804 -0.1643368 -0.3131850 -0.4902540 0.4063136 0.0000000 0.0000000
>
> # denominator of n1, n2, n3, diagonal of animals in LHS
> n_nu = a1 + a2 + a3
>
> # print
> n_nu
[1] 3.666667 4.000000 4.000000 4.666667 6.000000 6.000000 5.000000 5.000000
>
> # parents average fraction of breeding values
> pa = pa1 * (a1 / n_nu)
>
> # print
> pa
[1] 0.00000000 0.00000000 0.00000000 0.02812702 -0.01995143 0.02655816 -0.07775809 0.05431515
>
> # yield deviation fraction of breeding values
> yd = yd1 * (a2 / n_nu)
>
> # print
> yd
[,1]
[1,] 0.00000000
[2,] 0.00000000
[3,] 0.00000000
[4,] 0.03032093
[5,] -0.08407167
[6,] 0.08259500
[7,] -0.17170047
[8,] 0.12829953
>
> # progeny contribution
> pc1 = pc0 / a3
> pc1[is.nan(pc1) == TRUE] = 0
> pc1
[1] 0.21657807 -0.03754020 -0.08216841 -0.31318501 -0.49025399 0.40631358 0.00000000 0.00000000
>
> # progeny contribution fraction of breeding values
> pc = pc1 * (a3 / n_nu)
>
> # print
> pc
[1] 0.09844458 -0.01877010 -0.04108420 -0.06711107 -0.08170900 0.06771893 0.00000000 0.00000000
>
> # Progeny(Daughter) Yield Deviation(PYD, DYD)
>
> # n2 of progeny
> n2prog = a2 / (a1 + a2)
>
> # print
> n2prog
[1] 0.0000000 0.0000000 0.0000000 0.2727273 0.2000000 0.2000000 0.2000000 0.2000000
>
> # 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] 0.25343572 -0.01501608 0.02320667 -0.30625451 -0.34166831 0.26481591 0.00000000 0.00000000
> dyd_d
[1] 0.3818182 0.4000000 0.4000000 0.2000000 0.2000000 0.2000000 0.0000000 0.0000000
>
> # dyd
> dyd = dyd_n / dyd_d
> dyd[is.nan(dyd) == TRUE] = 0
>
> # print
> dyd
[1] 0.66376021 -0.03754020 0.05801667 -1.53127256 -1.70834154 1.32407954 0.00000000 0.00000000
>
> # breeding values and fractions
> stam_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
> stam_result
animal animal_bv rel acc sep pa yd pc sum_of_fr dyd
1 1 0.098444576 0.05781158 0.2404404 4.340941 0.00000000 0.00000000 0.09844458 0.098444576 0.66376021
2 2 -0.018770099 0.01580856 0.1257321 4.436646 0.00000000 0.00000000 -0.01877010 -0.018770099 -0.03754020
3 3 -0.041084203 0.08708243 0.2950973 4.272979 0.00000000 0.00000000 -0.04108420 -0.041084203 0.05801667
4 4 -0.008663123 0.14463969 0.3803153 4.136086 0.02812702 0.03032093 -0.06711107 -0.008663123 -1.53127256
5 5 -0.185732099 0.14378651 0.3791919 4.138148 -0.01995143 -0.08407167 -0.08170900 -0.185732099 -1.70834154
6 6 0.176872088 0.11543447 0.3397565 4.206104 0.02655816 0.08259500 0.06771893 0.176872088 1.32407954
7 7 -0.249458555 0.11628766 0.3410098 4.204075 -0.07775809 -0.17170047 0.00000000 -0.249458555 0.00000000
8 8 0.182614688 0.15527171 0.3940453 4.110300 0.05431515 0.12829953 0.00000000 0.182614688 0.00000000
>
> # 파일로 출력, 분리자는 ",", 따옴표 없고
> output_filename = gsub("[ -]", "", paste("stam_result_", Sys.Date(), ".csv"))
> write.table(stam_result, output_filename, sep=",", row.names = FALSE, quote = FALSE)
>
-
결과 파일 - stam_result_20200601.csv
animal,animal_bv,rel,acc,sep,pa,yd,pc,sum_of_fr,dyd
1,0.0984445757038791,0.0578115770821025,0.24044038155456,4.34094096462483,0,0,0.0984445757038788,0.0984445757038788,0.66376020776493
2,-0.0187700991008759,0.0158085581151151,0.125732088645322,4.43664612491212,0,0,-0.0187700991008722,-0.0187700991008722,-0.037540198201747
3,-0.0410842029270846,0.087082430924723,0.295097324496043,4.27297921613311,0,0,-0.0410842029270835,-0.0410842029270835,0.0580166699481492
4,-0.00866312266194402,0.144639692852924,0.38031525456248,4.13608584811069,0.0281270216296797,0.0303209293167937,-0.0671110736084169,-0.00866312266194338,-1.53127256021528
5,-0.18573209949465,0.143786506530158,0.379191912532635,4.13814812076572,-0.0199514340093202,-0.0840716676511125,-0.0817089978342196,-0.185732099494652,-1.70834153704798
6,0.176872087681303,0.11543446872744,0.33975648445238,4.20610397225879,0.0265581588676677,0.0825949990155542,0.0677189297980794,0.176872087681301,1.32407954321716
7,-0.249458554833631,0.116287655050207,0.341009757998517,4.20407503489125,-0.0777580888626375,-0.171700465970993,0,-0.24945855483363,0
8,0.182614687930696,0.155271707028943,0.39404531088308,4.11029997195109,0.0543151539016872,0.128299534029007,0,0.182614687930695,0
'Animal Breeding > R for Genetic Evaluation' 카테고리의 다른 글
| R을 이용한 반복 모형(Repeatability Model)의 육종가 구하기 프로그래밍 (0) | 2020.06.14 |
|---|---|
| R을 이용하여 반복 모형(repeatability model)에 대한 blupf90 분석 결과 정리하기 (0) | 2020.06.13 |
| R을 이용하여 단형질 개체 모형(single trait animal model)에 대한 blupf90 분석 결과 정리하기 (0) | 2020.06.03 |
| R을 이용한 혈통 파일 처리 (0) | 2020.04.10 |
| R로 단형질 임의회귀모형의 해 구하기 (1) | 2011.06.09 |