지금까지 여러 가지 모형의 육종가를 R로 구하였는데 모형이 복잡해지면서 R 프로그래밍 또한 복잡해 지고 있다. 각각의 모형에 알맞은 프로그램이 아니라 어느 정도 범용적으로 사용할 수 있는 프로그램을 만들어 보자. 본 내용은 Computational techniques in animal breeding, Ignacy Misztal (University of Georgia, Athens, USA)를 참고하였다. 위 책은 포트란 소스로 되어 있는데 그것을 R로 바꾸었다.
또한 농촌진흥청에 발행한 "동물 육종을 위한 육종가 추정 프로그램 작성"을 참고할 수 있다.
- single trait 2 fixed effect least square model
먼저 가장 기본적인 2개의 고정효과가 있는 최소제곱 모형을 다루어 보자.
자료
1 1 5
1 2 8
2 1 7
2 2 6
2 3 9
프로그램
# Single trait 2 fixed effects least squares model
# Date : 2020.08.05
# Reference : Computational techneques in animal breeding, Ignacy Misztal
# MASS packages
if (!("MASS" %in% installed.packages())) {
install.packages("MASS", dependencies = TRUE)
}
library(MASS)
# function : find address
find_addresses <- function() {
for (i in 1:neff) {
if (i == 1) {
address[i] <<- indata[i]
} else {
address[i] <<- sum(nlev[1:i - 1]) + indata[i]
}
}
}
# parameters
# number of effects
neff = 2
neff
# number of levels for each effect
nlev = c(2, 3)
nlev
# number of equations
neq = sum(nlev)
neq
# make a empty left hand side
lhs = matrix(rep(0, neq * neq), nrow = neq)
lhs
# make a empty right hand side
rhs = matrix(rep(0, neq))
rhs
# addresses
address = c(rep(0, neff))
address
# set working directory setwd(choose.dir())
setwd("D:/temp/00_prediction")
# print working directory
getwd()
# read data : effect1, effect2, y1
data = read.table("01_lsq_data.txt", header = FALSE, sep = "", col.names = c("e1",
"e2", "y1"), stringsAsFactors = FALSE)
data
# number of data
ndata = nrow(data)
ndata
# fill lhs and rhs
for (k in 1:ndata) {
indata = as.numeric(data[k, 1:neff])
y = as.numeric(data[k, neff + 1])
find_addresses()
for (i in 1:neff) {
for (j in 1:neff) {
lhs[address[i], address[j]] = lhs[address[i], address[j]] + 1
}
rhs[address[i]] = rhs[address[i]] + y
}
}
# print lhs, rhs
lhs
rhs
# generalised inverse of lhs
gi_lhs = ginv(lhs)
# print
gi_lhs
# solution
sol = gi_lhs %*% rhs
# print
sol
실행 결과
> # MASS packages
> if (!("MASS" %in% installed.packages())) {
+ install.packages("MASS", dependencies = TRUE)
+ }
> library(MASS)
> # function : find address
> find_addresses <- function() {
+
+ for (i in 1:neff) {
+ if (i == 1) {
+ address[i] <<- indata[i]
+ } else {
+ address[i] <<- sum(nlev[1:i - 1]) + indata[i]
+ }
+ }
+
+ }
> # number of effects
> neff = 2
> neff
[1] 2
> # number of levels for each effect
> nlev = c(2, 3)
> nlev
[1] 2 3
> # number of equations
> neq = sum(nlev)
> neq
[1] 5
> # make a empty left hand side
> lhs = matrix(rep(0, neq * neq), nrow = neq)
> lhs
[,1] [,2] [,3] [,4] [,5]
[1,] 0 0 0 0 0
[2,] 0 0 0 0 0
[3,] 0 0 0 0 0
[4,] 0 0 0 0 0
[5,] 0 0 0 0 0
> # make a empty right hand side
> rhs = matrix(rep(0, neq))
> rhs
[,1]
[1,] 0
[2,] 0
[3,] 0
[4,] 0
[5,] 0
> # addresses
> address = c(rep(0, neff))
> address
[1] 0 0
> # set working directory setwd(choose.dir())
> setwd("D:/temp/00_prediction")
> # print working directory
> getwd()
[1] "D:/temp/00_prediction"
> # read data : effect1, effect2, y1
> data = read.table("01_lsq_data.txt", header = FALSE, sep = "", col.names = c("e1",
+ "e2", "y1"), stringsAsFactors = FALSE)
> data
e1 e2 y1
1 1 1 5
2 1 2 8
3 2 1 7
4 2 2 6
5 2 3 9
> # number of data
> ndata = nrow(data)
> ndata
[1] 5
> # fill lhs and rhs
> for (k in 1:ndata) {
+ indata = as.numeric(data[k, 1:neff])
+ y = as.numeric(data[k, neff + 1])
+
+ find_addresses()
+
+ for (i in 1:neff) {
+ for (j in 1:neff) {
+ lhs[address[i], address[j]] = lhs[address[i], address[j]] + 1
+ }
+ rhs[address[i]] = rhs[address[i]] + y
+ }
+
+ }
> # print lhs, rhs
> lhs
[,1] [,2] [,3] [,4] [,5]
[1,] 2 0 1 1 0
[2,] 0 3 1 1 1
[3,] 1 1 2 0 0
[4,] 1 1 0 2 0
[5,] 0 1 0 0 1
> rhs
[,1]
[1,] 13
[2,] 22
[3,] 12
[4,] 14
[5,] 9
> # generalised inverse of lhs
> gi_lhs = ginv(lhs)
> # print
> gi_lhs
[,1] [,2] [,3] [,4] [,5]
[1,] 0.44 -0.16 -0.04 -0.04 0.36
[2,] -0.16 0.24 0.06 0.06 -0.04
[3,] -0.04 0.06 0.39 -0.11 -0.26
[4,] -0.04 0.06 -0.11 0.39 -0.26
[5,] 0.36 -0.04 -0.26 -0.26 0.84
> # solution
> sol = gi_lhs %*% rhs
> # print
> sol
[,1]
[1,] 4.4
[2,] 4.4
[3,] 1.6
[4,] 2.6
[5,] 4.6
>
-
single traits 1 fixed effect least squares model with covariable
다음은 공변이(covariable)가 있는 최소 제곱 모형이다.
자료
1 1 5
1 2 8
2 1 7
2 2 6
2 3 9
첫 열은 분류 변수, 둘째 열이 공변이 변수이다.
프로그램
# single traits 1 fixed effect least squares model with covariable
# Date : 2020.08.05
# Reference : Computational techneques in animal breeding, Ignacy Misztal
# MASS packages
if (!("MASS" %in% installed.packages())) {
install.packages("MASS", dependencies = TRUE)
}
library(MASS)
# function : find address
find_addresses <- function() {
for (i in 1:neff) {
if (efftype[i] == "effcross") {
weight_cov[i] <<- 1
if (i == 1) {
address[i] <<- indata[i]
} else {
address[i] <<- sum(nlev[1:i - 1]) + indata[i]
}
} else if (efftype[i] == "effcov") {
weight_cov[i] <<- indata[i]
if (nestedcov[i] == 0) {
if (i == 1) {
address[i] <<- 1
} else {
address[i] <<- sum(nlev[1:i - 1]) + 1
}
} else if (nestedcov[i] > 0 && nestedcov[i] < neff) {
if (i == 1) {
address[i] <<- indata[nestedcov[i]]
} else {
address[i] <<- sum(nlev[1:i - 1]) + indata[nestedcov[i]]
}
} else {
warning("wrong descriptin of nested covariable")
}
} else {
warning("unimplemented effect type")
}
}
}
# parameters
# number of effects
neff = 2
neff
# effect type
efftype = c("effcross", "effcov")
efftype
# nested covariable
nestedcov = c(0, 1)
nestedcov
# nlev
nlev = c(2, 2)
nlev
# set working directory setwd(choose.dir())
setwd("D:/temp/00_prediction")
# print working directory
getwd()
# read data : effect1, effect2, y1
data = read.table("02_lsqr_data.txt", header = FALSE, sep = "", col.names = c("eff1",
"eff2", "y1"), stringsAsFactors = FALSE)
data
# number of data
ndata = nrow(data)
ndata
# number of equation
neq = sum(nlev)
neq
# make a empty left hand side
lhs = matrix(rep(0, neq * neq), nrow = neq)
lhs
# make a empty right hand side
rhs = matrix(rep(0, neq))
rhs
# addresses
address = rep(0, neff)
address
# weight of covariable
weight_cov = rep(0, neff)
weight_cov
# fill lhs and rhs
for (k in 1:ndata) {
# k = 1
indata = as.numeric(data[k, 1:neff])
y = as.numeric(data[k, neff + 1])
find_addresses()
for (i in 1:neff) {
for (j in 1:neff) {
lhs[address[i], address[j]] = lhs[address[i], address[j]] + weight_cov[i] * weight_cov[j]
}
rhs[address[i]] = rhs[address[i]] + y * weight_cov[i]
}
}
# print lhs, rhs
lhs
rhs
# generalised inverse of lhs
gi_lhs = ginv(lhs)
# print
gi_lhs
# solution
sol = gi_lhs %*% rhs
# print
sol
실행 결과
> # MASS packages
> if (!("MASS" %in% installed.packages())) {
+ install.packages("MASS", dependencies = TRUE)
+ }
> library(MASS)
>
> # function : find address
> find_addresses <- function() {
+
+ for (i in 1:neff) {
+ if (efftype[i] == "effcross") {
+ weight_cov[i] <<- 1
+
+ if (i == 1) {
+ address[i] <<- indata[i]
+ } else {
+ address[i] <<- sum(nlev[1:i - 1]) + indata[i]
+ }
+
+
+ } else if (efftype[i] == "effcov") {
+ weight_cov[i] <<- indata[i]
+
+ if (nestedcov[i] == 0) {
+
+ if (i == 1) {
+ address[i] <<- indata[i]
+ } else {
+ address[i] <<- sum(nlev[1:i - 1]) + indata[i]
+ }
+ } else if (nestedcov[i] > 0 && nestedcov[i] < neff) {
+ if (i == 1) {
+ address[i] <<- indata[nestedcov[i]]
+ } else {
+ address[i] <<- sum(nlev[1:i - 1]) + indata[nestedcov[i]]
+ }
+
+ } else {
+ warning("wrong descriptin of nested covariable")
+ }
+
+ } else {
+ warning("unimplemented effect type")
+ }
+ }
+
+ }
>
> # parameters
>
> # number of effects
> neff = 2
> neff
[1] 2
>
> # effect type
> efftype = c("effcross", "effcov")
> efftype
[1] "effcross" "effcov"
>
> # nested covariable
> nestedcov = c(0, 1)
> nestedcov
[1] 0 1
>
> # nlev
> nlev = c(2, 2)
> nlev
[1] 2 2
>
> # set working directory setwd(choose.dir())
> setwd("D:/temp/00_prediction")
>
> # print working directory
> getwd()
[1] "D:/temp/00_prediction"
>
> # read data : effect1, effect2, y1
> data = read.table("02_lsqr_data.txt", header = FALSE, sep = "", col.names = c("eff1",
+ "eff2", "y1"), stringsAsFactors = FALSE)
> data
eff1 eff2 y1
1 1 1 5
2 1 2 8
3 2 1 7
4 2 2 6
5 2 3 9
>
> # number of data
> ndata = nrow(data)
> ndata
[1] 5
>
> # number of equation
> neq = sum(nlev)
> neq
[1] 4
>
> # make a empty left hand side
> lhs = matrix(rep(0, neq * neq), nrow = neq)
> lhs
[,1] [,2] [,3] [,4]
[1,] 0 0 0 0
[2,] 0 0 0 0
[3,] 0 0 0 0
[4,] 0 0 0 0
>
> # make a empty right hand side
> rhs = matrix(rep(0, neq))
> rhs
[,1]
[1,] 0
[2,] 0
[3,] 0
[4,] 0
>
> # addresses
> address = rep(0, neff)
> address
[1] 0 0
>
> # weight of covariable
> weight_cov = rep(0, neff)
> weight_cov
[1] 0 0
>
> # fill lhs and rhs
> for (k in 1:ndata) {
+ # k = 1
+ indata = as.numeric(data[k, 1:neff])
+ y = as.numeric(data[k, neff + 1])
+
+ find_addresses()
+
+ for (i in 1:neff) {
+ for (j in 1:neff) {
+ lhs[address[i], address[j]] = lhs[address[i], address[j]] + weight_cov[i] * weight_cov[j]
+ }
+ rhs[address[i]] = rhs[address[i]] + y * weight_cov[i]
+ }
+
+ }
>
> # print lhs, rhs
> lhs
[,1] [,2] [,3] [,4]
[1,] 2 0 3 0
[2,] 0 3 0 6
[3,] 3 0 5 0
[4,] 0 6 0 14
> rhs
[,1]
[1,] 13
[2,] 22
[3,] 21
[4,] 46
>
> # generalised inverse of lhs
> gi_lhs = ginv(lhs)
>
> # print
> gi_lhs
[,1] [,2] [,3] [,4]
[1,] 5 0.000000 -3 0.0
[2,] 0 2.333333 0 -1.0
[3,] -3 0.000000 2 0.0
[4,] 0 -1.000000 0 0.5
>
> # solution
> sol = gi_lhs %*% rhs
>
> # print
> sol
[,1]
[1,] 2.000000
[2,] 5.333333
[3,] 3.000000
[4,] 1.000000
>
-
Inverse of Numerator Relationship Matirx
animal model을 이용하기 위해선 개체들 사이의 혈연 계수 행렬 또는 혈연 계수 행렬의 역행렬을 구해야 한다. 다음은 혈연 계수 행렬의 역행렬을 구하는 프로그램이다.
자료
1 0 0
2 0 0
3 0 0
4 1 2
5 3 2
6 1 5
7 3 4
8 1 7
프로그램
# inverse of numerator relationship matrix
# Date : 2020.08.05
# Reference : Computational techneques in animal breeding, Ignacy Misztal
# function to make inverse of numerator relationship matrix
ainv = function() {
w = c(1, -0.5, -0.5)
res = c(2, 4/3, 1, 0)
for (di in 1:nped) {
p = as.numeric(pedi[di,])
i = 1
if(p[2] == 0){ i = i + 1 }
if(p[3] == 0){ i = i + 1 }
mendel = res[i]
for (k in 1:3) {
for (l in 1:3) {
if (p[k] != 0 && p[l] != 0) {
a[p[k], p[l]] <<- a[p[k], p[l]] + w[k] * w[l] * mendel
}
}
}
}
}
# set working directory setwd(choose.dir())
setwd("D:/temp/00_prediction")
# print working directory
getwd()
# read pedigree data
pedi = read.table("03_ainv_data.txt", header = FALSE, sep = "", col.names = c("a",
"s", "d"), stringsAsFactors = FALSE)
pedi
# number of pedi
nped = nrow(pedi)
nped
# inverse matrix of NRM(Nemerator Relationship Matrix)
a = matrix(c(0), nrow = nped, ncol = nped)
a
# run function ainv
ainv()
# print NRM
a
실형 결과
> # function to make inverse of numerator relationship matrix
> ainv = function() {
+
+ w = c(1, -0.5, -0.5)
+ res = c(2, 4/3, 1, 0)
+
+ for (di in 1:nped) {
+
+ p = as.numeric(pedi[di,])
+
+ i = 1
+ if(p[2] == 0){ i = i + 1 }
+ if(p[3] == 0){ i = i + 1 }
+
+ mendel = res[i]
+
+ for (k in 1:3) {
+ for (l in 1:3) {
+ if (p[k] != 0 && p[l] != 0) {
+ a[p[k], p[l]] <<- a[p[k], p[l]] + w[k] * w[l] * mendel
+ }
+ }
+ }
+
+ }
+
+ }
>
> # set working directory setwd(choose.dir())
> setwd("D:/temp/00_prediction")
>
> # print working directory
> getwd()
[1] "D:/temp/00_prediction"
>
> # read pedigree data
> pedi = read.table("03_ainv_data.txt", header = FALSE, sep = "", col.names = c("a",
+ "s", "d"), stringsAsFactors = FALSE)
> pedi
a s d
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
>
> # number of pedi
> nped = nrow(pedi)
> nped
[1] 8
>
> # inverse matrix of NRM(Nemerator Relationship Matrix)
> a = matrix(c(0), nrow = nped, ncol = nped)
> a
[,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 0 0 0 0 0
[5,] 0 0 0 0 0 0 0 0
[6,] 0 0 0 0 0 0 0 0
[7,] 0 0 0 0 0 0 0 0
[8,] 0 0 0 0 0 0 0 0
>
> # run function ainv
> ainv()
>
> # print NRM
> a
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 2.5 0.5 0.0 -1.0 0.5 -1 0.5 -1
[2,] 0.5 2.0 0.5 -1.0 -1.0 0 0.0 0
[3,] 0.0 0.5 2.0 0.5 -1.0 0 -1.0 0
[4,] -1.0 -1.0 0.5 2.5 0.0 0 -1.0 0
[5,] 0.5 -1.0 -1.0 0.0 2.5 -1 0.0 0
[6,] -1.0 0.0 0.0 0.0 -1.0 2 0.0 0
[7,] 0.5 0.0 -1.0 -1.0 0.0 0 2.5 -1
[8,] -1.0 0.0 0.0 0.0 0.0 0 -1.0 2
>
-
single trait animal model
다음은 단형질 개체 모형이다.
표현형 자료
1 2 8
1 1 5
2 1 7
2 2 6
2 3 9
1열: 분류변수, 2열 : 개체, 3열 : 표현형 자료
혈통 자료
5 0 0
1 3 4
2 3 4
3 5 0
4 0 0
프로그램
# single trait animal model
# Date : 2020.08.05
# Reference : Computational techneques in animal breeding, Ignacy Misztal
# MASS packages
if (!("MASS" %in% installed.packages())) {
install.packages("MASS", dependencies = TRUE)
}
library(MASS)
# function : find address
find_addresses <- function() {
for (i in 1:neff) {
if (i == 1) {
address[i] <<- indata[i]
} else {
address[i] <<- sum(nlev[1:i - 1]) + indata[i]
}
}
}
# address for given level l of effect e
address1 <- function(e, l){
ifelse(e == 1, l, sum(nlev[1 : e - 1]) + l)
return(ifelse(e == 1, l, sum(nlev[1 : e - 1]) + l))
}
# function to make inverse of numerator relationship matrix
ainv = function(eff) {
w = c(1, -0.5, -0.5)
res = c(2, 4/3, 1, 0)
for (di in 1:nped) {
p = as.numeric(pedi[di,])
i = 1
if(p[2] == 0){ i = i + 1 }
if(p[3] == 0){ i = i + 1 }
mendel = res[i]
for (k in 1:3) {
for (l in 1:3) {
if (p[k] != 0 && p[l] != 0) {
m = address1(eff, p[k])
n = address1(eff, p[l])
a[p[k], p[l]] <<- a[p[k], p[l]] + w[k] * w[l] * mendel
lhs[m, n] <<- lhs[m, n] + w[k] * w[l] * mendel / var_a
}
}
}
}
}
# when random effect is diagonal
diag <- function(eff){
for(i in 1 : nlev[eff]){
m = address1(eff, i)
lhs[m, m] <<- lhs[m, m] + 1 / var_a
}
}
# parameters
# number of effects
neff = 2
neff
# nlev
nlev = c(2, 5)
nlev
# random effect number
addeff = 2
addeff
# type of random effect
random = "diag"
random
# additive genetic variance
var_a = 2
var_a
# residual variance
var_e = 10
var_e
# number of equation
neq = sum(nlev)
neq
# empty left hand side
lhs = matrix(c(0), nrow = neq, ncol=neq)
lhs
# empty right hand side
rhs = matrix(c(0), nrow = neq)
rhs
# addresses
address = c(rep(0, neff))
address
# set working directory setwd(choose.dir())
setwd("D:/temp/00_prediction")
# print working directory
getwd()
# read data : effect1, effect2, y1
data = read.table("04_stam_data.txt", header = FALSE, sep = "", col.names = c("e1",
"e2", "y1"), stringsAsFactors = FALSE)
data
# number of data
ndata = nrow(data)
ndata
# read pedigree : animal, sire, dam
pedi = read.table("04_stam_pedi.txt", header = FALSE, sep = "", col.names = c("animal",
"sire", "dam"), stringsAsFactors = FALSE)
pedi
# number of pedigree
nped = nrow(pedi)
nped
# inverse matrix of NRM(Nemerator Relationship Matrix)
a = matrix(c(0), nrow = nped, ncol = nped)
a
# fill lhs and rhs
for (k in 1 : ndata) {
# k = 1
indata = as.numeric(data[k, 1:neff])
y = as.numeric(data[k, neff + 1])
find_addresses()
for (i in 1:neff) {
for (j in 1:neff) {
lhs[address[i], address[j]] = lhs[address[i], address[j]] + 1/var_e
}
rhs[address[i]] = rhs[address[i]] + y/var_e
}
}
# print lhs, rhs
lhs
rhs
# process random effect
for (i in 1:neff) {
if (i == addeff) {
if (random == "add") {
# add the inverse of NRM to lhs
ainv(i)
} else if (random == "diag") {
# add the inverse of var_a to the diagonal of lhs
diag(i)
} else {
warning("Choose add or diag as a type of random effect")
}
}
}
# print
lhs
rhs
# generalised inverse of lhs
gi_lhs = ginv(lhs)
# print
gi_lhs
# solution
sol = gi_lhs %*% rhs
# print
sol
실행 결과
> # MASS packages
> if (!("MASS" %in% installed.packages())) {
+ install.packages("MASS", dependencies = TRUE)
+ }
> library(MASS)
>
> # function : find address
> find_addresses <- function() {
+
+ for (i in 1:neff) {
+ if (i == 1) {
+ address[i] <<- indata[i]
+ } else {
+ address[i] <<- sum(nlev[1:i - 1]) + indata[i]
+ }
+ }
+
+ }
>
> # address for given level l of effect e
> address1 <- function(e, l){
+ ifelse(e == 1, l, sum(nlev[1 : e - 1]) + l)
+ return(ifelse(e == 1, l, sum(nlev[1 : e - 1]) + l))
+ }
>
> # function to make inverse of numerator relationship matrix
> ainv = function(eff) {
+
+ w = c(1, -0.5, -0.5)
+ res = c(2, 4/3, 1, 0)
+
+ for (di in 1:nped) {
+
+ p = as.numeric(pedi[di,])
+
+ i = 1
+ if(p[2] == 0){ i = i + 1 }
+ if(p[3] == 0){ i = i + 1 }
+
+ mendel = res[i]
+
+ for (k in 1:3) {
+ for (l in 1:3) {
+ if (p[k] != 0 && p[l] != 0) {
+ m = address1(eff, p[k])
+ n = address1(eff, p[l])
+ a[p[k], p[l]] <<- a[p[k], p[l]] + w[k] * w[l] * mendel
+ lhs[m, n] <<- lhs[m, n] + w[k] * w[l] * mendel / var_a
+ }
+ }
+ }
+
+ }
+
+ }
>
> # when random effect is diagonal
> diag <- function(eff){
+ for(i in 1 : nlev[eff]){
+ m = address1(eff, i)
+ lhs[m, m] <<- lhs[m, m] + 1 / var_a
+ }
+ }
>
> # parameters
>
> # number of effects
> neff = 2
> neff
[1] 2
>
> # nlev
> nlev = c(2, 5)
> nlev
[1] 2 5
>
> # random effect number
> addeff = 2
> addeff
[1] 2
>
> # type of random effect
> random = "diag"
> random
[1] "diag"
>
> # additive genetic variance
> var_a = 2
> var_a
[1] 2
>
> # residual variance
> var_e = 10
> var_e
[1] 10
>
> # number of equation
> neq = sum(nlev)
> neq
[1] 7
>
> # empty left hand side
> lhs = matrix(c(0), nrow = neq, ncol=neq)
> lhs
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 0 0 0 0 0 0 0
[2,] 0 0 0 0 0 0 0
[3,] 0 0 0 0 0 0 0
[4,] 0 0 0 0 0 0 0
[5,] 0 0 0 0 0 0 0
[6,] 0 0 0 0 0 0 0
[7,] 0 0 0 0 0 0 0
>
> # empty right hand side
> rhs = matrix(c(0), nrow = neq)
> rhs
[,1]
[1,] 0
[2,] 0
[3,] 0
[4,] 0
[5,] 0
[6,] 0
[7,] 0
>
> # addresses
> address = c(rep(0, neff))
> address
[1] 0 0
>
> # set working directory setwd(choose.dir())
> setwd("D:/temp/00_prediction")
>
> # print working directory
> getwd()
[1] "D:/temp/00_prediction"
>
> # read data : effect1, effect2, y1
> data = read.table("04_stam_data.txt", header = FALSE, sep = "", col.names = c("e1",
+ "e2", "y1"), stringsAsFactors = FALSE)
> data
e1 e2 y1
1 1 2 8
2 1 1 5
3 2 1 7
4 2 2 6
5 2 3 9
>
> # number of data
> ndata = nrow(data)
> ndata
[1] 5
>
> # read pedigree : animal, sire, dam
> pedi = read.table("04_stam_pedi.txt", header = FALSE, sep = "", col.names = c("animal",
+ "sire", "dam"), stringsAsFactors = FALSE)
> pedi
animal sire dam
1 5 0 0
2 1 3 4
3 2 3 4
4 3 5 0
5 4 0 0
>
> # number of pedigree
> nped = nrow(pedi)
> nped
[1] 5
>
> # inverse matrix of NRM(Nemerator Relationship Matrix)
> a = matrix(c(0), nrow = nped, ncol = nped)
> a
[,1] [,2] [,3] [,4] [,5]
[1,] 0 0 0 0 0
[2,] 0 0 0 0 0
[3,] 0 0 0 0 0
[4,] 0 0 0 0 0
[5,] 0 0 0 0 0
>
> # fill lhs and rhs
>
> for (k in 1 : ndata) {
+ # k = 1
+ indata = as.numeric(data[k, 1:neff])
+ y = as.numeric(data[k, neff + 1])
+
+ find_addresses()
+
+ for (i in 1:neff) {
+ for (j in 1:neff) {
+ lhs[address[i], address[j]] = lhs[address[i], address[j]] + 1/var_e
+ }
+ rhs[address[i]] = rhs[address[i]] + y/var_e
+ }
+
+ }
>
> # print lhs, rhs
> lhs
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 0.2 0.0 0.1 0.1 0.0 0 0
[2,] 0.0 0.3 0.1 0.1 0.1 0 0
[3,] 0.1 0.1 0.2 0.0 0.0 0 0
[4,] 0.1 0.1 0.0 0.2 0.0 0 0
[5,] 0.0 0.1 0.0 0.0 0.1 0 0
[6,] 0.0 0.0 0.0 0.0 0.0 0 0
[7,] 0.0 0.0 0.0 0.0 0.0 0 0
> rhs
[,1]
[1,] 1.3
[2,] 2.2
[3,] 1.2
[4,] 1.4
[5,] 0.9
[6,] 0.0
[7,] 0.0
>
> # process random effect
> for (i in 1:neff) {
+ if (i == addeff) {
+
+ if (random == "add") {
+
+ # add the inverse of NRM to lhs
+ ainv(i)
+
+ } else if (random == "diag") {
+
+ # add the inverse of var_a to the diagonal of lhs
+ diag(i)
+
+ } else {
+ warning("Choose add or diag as a type of random effect")
+ }
+ }
+ }
>
> # print
> lhs
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 0.2 0.0 0.1 0.1 0.0 0.0 0.0
[2,] 0.0 0.3 0.1 0.1 0.1 0.0 0.0
[3,] 0.1 0.1 0.7 0.0 0.0 0.0 0.0
[4,] 0.1 0.1 0.0 0.7 0.0 0.0 0.0
[5,] 0.0 0.1 0.0 0.0 0.6 0.0 0.0
[6,] 0.0 0.0 0.0 0.0 0.0 0.5 0.0
[7,] 0.0 0.0 0.0 0.0 0.0 0.0 0.5
> rhs
[,1]
[1,] 1.3
[2,] 2.2
[3,] 1.2
[4,] 1.4
[5,] 0.9
[6,] 0.0
[7,] 0.0
>
> # generalised inverse of lhs
> gi_lhs = ginv(lhs)
>
> # print
> gi_lhs
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 5.9444444 0.6666667 -0.9444444 -0.9444444 -0.1111111 0 0
[2,] 0.6666667 4.0000000 -0.6666667 -0.6666667 -0.6666667 0 0
[3,] -0.9444444 -0.6666667 1.6587302 0.2301587 0.1111111 0 0
[4,] -0.9444444 -0.6666667 0.2301587 1.6587302 0.1111111 0 0
[5,] -0.1111111 -0.6666667 0.1111111 0.1111111 1.7777778 0 0
[6,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 2 0
[7,] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0 2
>
> # solution
> sol = gi_lhs %*% rhs
>
> # print
> sol
[,1]
[1,] 6.638888889
[2,] 7.333333333
[3,] -0.281746032
[4,] 0.003968254
[5,] 0.277777778
[6,] 0.000000000
[7,] 0.000000000
>
-
multiple traits least squares model
다형질 개체 모형을 다루어야 하는데 우선 최소 제곱 모형을 다루고 후에 다형질 개체 모형을 다루도록 하자
자료
2 3 2 3 9 1
2 2 2 2 6 5
2 1 2 1 7 3
1 2 1 2 8 4
1 1 1 1 5 2
1, 2열 : 첫째 형질의 2개의 fixed effects
3, 4열 : 둘째 형질의 2개의 fixed effects
5, 6열 : 첫째, 둘째 형질의 표현형
프로그램
# multiple traits least squares model
# Date : 2020.08.06.
# Reference : Computational techneques in animal breeding, Ignacy Misztal
# MASS packages
if (!("MASS" %in% installed.packages())) {
install.packages("MASS", dependencies = TRUE)
}
library(MASS)
# return data in indata
datum <- function(ef, tr) {
return(indata[ef + (tr - 1) * neff])
}
find_rinv <- function() {
tempr = r
for (i in 1:ntrait) {
if (y[i] == 0) {
tempr[i, ] = 0
tempr[, i] = 0
}
}
rinv <<- ginv(tempr)
}
# function : find address
find_addresses <- function() {
for (i in 1:neff) {
for (j in 1:ntrait) {
# i=2 j=1
if (datum(i, j) == miss) {
address[i, j] <<- 0
weight_cov[i, j] <<- 0
next
}
baseaddr = ifelse(i == 1, j, sum(nlev[1:i - 1]) * ntrait + j)
if (effecttype[i] == "effcross") {
weight_cov[i, j] <<- 1
address[i, j] <<- baseaddr + (datum(i, j) - 1) * ntrait
} else if (effecttype[i] == "effcov") {
weight_cov[i, j] <<- datum(i, j)
if (nestedcov[i] == 0) {
address[i, j] <<- baseaddr
} else if (nestedcov[i] > 0 && nestedcov[i] < neff) {
address[i, j] <<- baseaddr + (datum(nestedcov[i], j) - 1) * ntrait
} else {
warning("wrong description of nested covariable")
}
} else {
warning("unimplemented effect")
}
}
}
}
# parameters
# number of traits
ntrait = 2
ntrait
# number of effects
neff = 2
neff
# type of effects
effecttype = c("effcross", "effcross")
effecttype
# whether covariable effect is nested or not
nestedcov = c(0, 0)
nestedcov
# number of levels
nlev = c(2, 3)
nlev
# missing values
miss = 0
miss
# residual variace
r = matrix(c(1, 2, 2, 5), nrow = 2)
r
# empty inverse of r
rinv = matrix(c(0), nrow = ntrait, ncol = ntrait)
rinv
# number of equation
neq = ntrait * sum(nlev)
neq
# make empty lhs and rhs
lhs = matrix(c(0), nrow = neq, ncol = neq)
lhs
rhs = matrix(c(0), nrow = neq)
rhs
# addresses
address = matrix(c(0), nrow=neff, ncol=ntrait)
address
# weight of covariable
weight_cov = matrix(c(0), nrow=neff, ncol=ntrait)
weight_cov
# set working directory setwd(choose.dir())
setwd("D:/temp/00_prediction")
# print working directory
getwd()
# read data : t1e1, t1e2, t2e1, t2e2, y1, y2
data = read.table("06_mt_lsq_data.txt", header = FALSE, sep = "", col.names = c("t1e1", "t1e2", "t2e1", "t2e2", "y1", "y2"), stringsAsFactors = FALSE)
data
# number of data
ndata = nrow(data)
ndata
# loop data
for (di in 1:ndata) {
# di = 5
indata = as.numeric(data[di, 1:(ntrait * neff)])
y = as.numeric(data[di, (ntrait * neff + 1):(ntrait * (neff + 1))])
# address and weight_cov
find_addresses()
# inverse of r
find_rinv()
for (i in 1:neff) {
for (j in 1:neff) {
for (k in 1:ntrait) {
for (l in 1:ntrait) {
lhs[address[i, k], address[j, l]] = lhs[address[i, k], address[j, l]] + weight_cov[i, k] * weight_cov[j, l] * rinv[k, l]
}
}
}
for (k in 1:ntrait) {
for (l in 1:ntrait) {
rhs[address[i, k]] = rhs[address[i, k]] + weight_cov[i, k] * rinv[k, l] * y[l]
}
}
}
}
# print
lhs
rhs
# generalised inverse of lhs
gi_lhs = ginv(lhs)
# print
gi_lhs
# solution
sol = gi_lhs %*% rhs
# print
sol
실형 결과
> # MASS packages
> if (!("MASS" %in% installed.packages())) {
+ install.packages("MASS", dependencies = TRUE)
+ }
> library(MASS)
>
> # return data in indata
> datum <- function(ef, tr) {
+ return(indata[ef + (tr - 1) * neff])
+ }
>
> find_rinv <- function() {
+ tempr = r
+ for (i in 1:ntrait) {
+ if (y[i] == 0) {
+ tempr[i, ] = 0
+ tempr[, i] = 0
+ }
+ }
+ rinv <<- ginv(tempr)
+ }
>
> # function : find address
> find_addresses <- function() {
+
+ for (i in 1:neff) {
+ for (j in 1:ntrait) {
+
+ # i=2 j=1
+
+ if (datum(i, j) == miss) {
+ address[i, j] <<- 0
+ weight_cov[i, j] <<- 0
+ next
+ }
+
+ baseaddr = ifelse(i == 1, j, sum(nlev[1:i - 1]) * ntrait + j)
+
+ if (effecttype[i] == "effcross") {
+ weight_cov[i, j] <<- 1
+ address[i, j] <<- baseaddr + (datum(i, j) - 1) * ntrait
+ } else if (effecttype[i] == "effcov") {
+ weight_cov[i, j] <<- datum(i, j)
+ if (nestedcov[i] == 0) {
+ address[i, j] <<- baseaddr
+ } else if (nestedcov[i] > 0 && nestedcov[i] < neff) {
+ address[i, j] <<- baseaddr + (datum(nestedcov[i], j) - 1) * ntrait
+ } else {
+ warning("wrong description of nested covariable")
+ }
+ } else {
+ warning("unimplemented effect")
+ }
+ }
+ }
+
+ }
>
>
> # parameters
>
> # number of traits
> ntrait = 2
> ntrait
[1] 2
>
> # number of effects
> neff = 2
> neff
[1] 2
>
> # type of effects
> effecttype = c("effcross", "effcross")
> effecttype
[1] "effcross" "effcross"
>
> # whether covariable effect is nested or not
> nestedcov = c(0, 0)
> nestedcov
[1] 0 0
>
> # number of levels
> nlev = c(2, 3)
> nlev
[1] 2 3
>
> # missing values
> miss = 0
> miss
[1] 0
>
> # residual variace
> r = matrix(c(1, 2, 2, 5), nrow = 2)
> r
[,1] [,2]
[1,] 1 2
[2,] 2 5
>
> # empty inverse of r
> rinv = matrix(c(0), nrow = ntrait, ncol = ntrait)
> rinv
[,1] [,2]
[1,] 0 0
[2,] 0 0
>
> # number of equation
> neq = ntrait * sum(nlev)
> neq
[1] 10
>
> # make empty lhs and rhs
> lhs = matrix(c(0), nrow = neq, ncol = neq)
> lhs
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0 0 0 0 0 0 0 0 0 0
[2,] 0 0 0 0 0 0 0 0 0 0
[3,] 0 0 0 0 0 0 0 0 0 0
[4,] 0 0 0 0 0 0 0 0 0 0
[5,] 0 0 0 0 0 0 0 0 0 0
[6,] 0 0 0 0 0 0 0 0 0 0
[7,] 0 0 0 0 0 0 0 0 0 0
[8,] 0 0 0 0 0 0 0 0 0 0
[9,] 0 0 0 0 0 0 0 0 0 0
[10,] 0 0 0 0 0 0 0 0 0 0
>
> rhs = matrix(c(0), nrow = neq)
> rhs
[,1]
[1,] 0
[2,] 0
[3,] 0
[4,] 0
[5,] 0
[6,] 0
[7,] 0
[8,] 0
[9,] 0
[10,] 0
>
> # addresses
> address = matrix(c(0), nrow=neff, ncol=ntrait)
> address
[,1] [,2]
[1,] 0 0
[2,] 0 0
>
> # weight of covariable
> weight_cov = matrix(c(0), nrow=neff, ncol=ntrait)
> weight_cov
[,1] [,2]
[1,] 0 0
[2,] 0 0
>
> # set working directory setwd(choose.dir())
> setwd("D:/temp/00_prediction")
>
> # print working directory
> getwd()
[1] "D:/temp/00_prediction"
>
> # read data : t1e1, t1e2, t2e1, t2e2, y1, y2
> data = read.table("06_mt_lsq_data.txt", header = FALSE, sep = "", col.names = c("t1e1", "t1e2", "t2e1", "t2e2", "y1", "y2"), stringsAsFactors = FALSE)
> data
t1e1 t1e2 t2e1 t2e2 y1 y2
1 2 3 2 3 9 1
2 2 2 2 2 6 5
3 2 1 2 1 7 3
4 1 2 1 2 8 4
5 1 1 1 1 5 2
>
> # number of data
> ndata = nrow(data)
> ndata
[1] 5
>
>
> # loop data
> for (di in 1:ndata) {
+ # di = 5
+ indata = as.numeric(data[di, 1:(ntrait * neff)])
+ y = as.numeric(data[di, (ntrait * neff + 1):(ntrait * (neff + 1))])
+
+ # address and weight_cov
+ find_addresses()
+
+ # inverse of r
+ find_rinv()
+
+ for (i in 1:neff) {
+ for (j in 1:neff) {
+ for (k in 1:ntrait) {
+ for (l in 1:ntrait) {
+ lhs[address[i, k], address[j, l]] = lhs[address[i, k], address[j, l]] + weight_cov[i, k] * weight_cov[j, l] * rinv[k, l]
+ }
+ }
+ }
+ for (k in 1:ntrait) {
+ for (l in 1:ntrait) {
+ rhs[address[i, k]] = rhs[address[i, k]] + weight_cov[i, k] * rinv[k, l] * y[l]
+ }
+ }
+ }
+ }
>
> # print
> lhs
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 10 -4 0 0 5 -2 5 -2 0 0
[2,] -4 2 0 0 -2 1 -2 1 0 0
[3,] 0 0 15 -6 5 -2 5 -2 5 -2
[4,] 0 0 -6 3 -2 1 -2 1 -2 1
[5,] 5 -2 5 -2 10 -4 0 0 0 0
[6,] -2 1 -2 1 -4 2 0 0 0 0
[7,] 5 -2 5 -2 0 0 10 -4 0 0
[8,] -2 1 -2 1 0 0 -4 2 0 0
[9,] 0 0 5 -2 0 0 0 0 5 -2
[10,] 0 0 -2 1 0 0 0 0 -2 1
> rhs
[,1]
[1,] 53
[2,] -20
[3,] 92
[4,] -35
[5,] 50
[6,] -19
[7,] 52
[8,] -19
[9,] 43
[10,] -17
>
> # generalised inverse of lhs
> gi_lhs = ginv(lhs)
>
> # print
> gi_lhs
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0.44 0.88 -0.16 -0.32 -0.04 -0.08 -0.04 -0.08 0.36 0.72
[2,] 0.88 2.20 -0.32 -0.80 -0.08 -0.20 -0.08 -0.20 0.72 1.80
[3,] -0.16 -0.32 0.24 0.48 0.06 0.12 0.06 0.12 -0.04 -0.08
[4,] -0.32 -0.80 0.48 1.20 0.12 0.30 0.12 0.30 -0.08 -0.20
[5,] -0.04 -0.08 0.06 0.12 0.39 0.78 -0.11 -0.22 -0.26 -0.52
[6,] -0.08 -0.20 0.12 0.30 0.78 1.95 -0.22 -0.55 -0.52 -1.30
[7,] -0.04 -0.08 0.06 0.12 -0.11 -0.22 0.39 0.78 -0.26 -0.52
[8,] -0.08 -0.20 0.12 0.30 -0.22 -0.55 0.78 1.95 -0.52 -1.30
[9,] 0.36 0.72 -0.04 -0.08 -0.26 -0.52 -0.26 -0.52 0.84 1.68
[10,] 0.72 1.80 -0.08 -0.20 -0.52 -1.30 -0.52 -1.30 1.68 4.20
>
> # solution
> sol = gi_lhs %*% rhs
>
> # print
> sol
[,1]
[1,] 4.4
[2,] 1.0
[3,] 4.4
[4,] 2.0
[5,] 1.6
[6,] 1.0
[7,] 2.6
[8,] 3.0
[9,] 4.6
[10,] -1.0
>
-
multiple trait animal model
마지막으로 다형질 개체 모형이다. 이 프로그램은 missing record, unequal desing matrices 까지 처리 가능한다.
자료
2 3 2 3 9 1
2 2 2 2 6 5
2 1 2 1 7 3
1 2 1 2 8 4
1 1 1 1 5 2
1, 2열 : 첫째 형질의 2개의 fixed effects
3, 4열 : 둘째 형질의 2개의 fixed effects
5, 6열 : 첫째, 둘째 형질의 표현형
혈통
1 3 4
2 3 4
3 5 0
4 0 0
5 0 0
프로그램
# multiple trait animal model
# Date : 2020.08.06.
# Reference : Computational techneques in animal breeding, Ignacy Misztal
# MASS packages
if (!("MASS" %in% installed.packages())) {
install.packages("MASS", dependencies = TRUE)
}
library(MASS)
# function : find address
find_addresses <- function() {
for (i in 1:neff) {
for (j in 1:ntrait) {
# i=2 j=1
if (datum(i, j) == miss) {
address[i, j] <<- 0
weight_cov[i, j] <<- 0
next
}
baseaddr = ifelse(i == 1, j, sum(nlev[1:i - 1]) * ntrait + j)
if (effecttype[i] == "effcross") {
weight_cov[i, j] <<- 1
address[i, j] <<- baseaddr + (datum(i, j) - 1) * ntrait
} else if (effecttype[i] == "effcov") {
weight_cov[i, j] <<- datum(i, j)
if (nestedcov[i] == 0) {
address[i, j] <<- baseaddr
} else if (nestedcov[i] > 0 && nestedcov[i] < neff) {
address[i, j] <<- baseaddr + (datum(nestedcov[i], j) - 1) * ntrait
} else {
warning("wrong description of nested covariable")
}
} else {
warning("unimplemented effect")
}
}
}
}
# return data in indata
datum <- function(ef, tr) {
return(indata[ef + (tr - 1) * neff])
}
# inverse of residual variance and covariance matrix
find_rinv <- function() {
tempr = r
for (i in 1:ntrait) {
if (y[i] == miss) {
tempr[i, ] = 0
tempr[, i] = 0
}
}
rinv <<- ginv(tempr)
}
# address for given level l of effect e and trait t
address1 <- function(e, l, t){
return(ifelse(e == 1, (l - 1) * ntrait + t, sum(nlev[1 : e - 1]) * ntrait + (l - 1) * ntrait + t))
}
# setup g
setup_g <- function(){
for(i in 1 : neff){
if (randomnumb[i] != 0){
g[,,i] <<- ginv(g0[,,i])
}
}
}
# add a genetic variance to the diagonal of left hand side
add_g_diag <- function(eff){
for(i in 1 : nlev[eff]){
for(j in 0 : (randomnumb[eff] - 1)){
for(k in 0 : (randomnumb[eff] - 1)){
for(t1 in 1 : ntrait){
for(t2 in 1 : ntrait){
m = address1(eff + j, i, t1)
l = address1(eff + k, i, t2)
#print(paste("i=",i,"j=",j,"k=",k,"t1=",t1,"t2=",t2,"m=",m,"l=",l))
lhs[m, l] <<- lhs[m, l] + g[t1 + j * ntrait, t2 + k * ntrait, eff]
}
}
}
}
}
}
# add genetic variacne matrix * NRN to left hand side
add_g_add <- function(type, eff){
if(type == 'g_A'){
w = c(1, -0.5, -0.5)
res = c(2, 4/3, 1, 0)
} else if(type == 'g_As'){
w = c(1, -0.5, -0.25)
res = c(16/11, 4/3, 16/15, 1)
}
for (di in 1:nped) {
p = as.numeric(pedi[di,])
ri = 1
if(p[2] == 0){ ri = ri + 1 }
if(p[3] == 0){ ri = ri + 1 }
mendel = res[ri]
for(i in 0 : (randomnumb[eff] - 1)){
for(j in 0 : (randomnumb[eff] - 1)){
for(t1 in 1 : ntrait){
for(t2 in 1 : ntrait){
for (k in 1 : 3) {
for (l in 1 : 3) {
if (p[k] != 0 && p[l] != 0) {
m = address1(eff + i, p[k], t1)
n = address1(eff + j, p[l], t2)
a[p[k], p[l]] <<- a[p[k], p[l]] + w[k] * w[l] * mendel
lhs[m, n] <<- lhs[m, n] + g[t1 + i * ntrait, t2 + j * ntrait, eff] * w[k] * w[l] * mendel
}
}
}
}
}
}
}
}
}
# parameters
# number of traits
ntrait = 2
ntrait
# number of effects
neff = 2
neff
# type of effects
effecttype = c("effcross", "effcross")
effecttype
# whether covariable effect is nested or not
nestedcov = c(0, 0)
nestedcov
# number of levels
nlev = c(2, 5)
nlev
# missing values
miss = 0
miss
# type of random effect : g_fixed, g_diag, g_A, g_As
randomtype = c("g_fixed", "g_A")
randomtype
# number of correlated effects per effect
randomnumb = c(0, 1)
randomnumb
# genetic varinace and covariance matrix
g0 = array(c(0), dim =c(ntrait, ntrait, neff))
g0[,,2] = matrix(c(2, -1, -1, 1), nrow = ntrait)
g0
# matrix for inverse of genetic variance matrix
g = array(c(0), dim =c(ntrait, ntrait, neff))
g
setup_g()
# residual variace
r = matrix(c(1, 0, 0, 1), nrow = 2)
r
# empty inverse of r
rinv = matrix(c(0), nrow = ntrait, ncol = ntrait)
rinv
# number of equations
neq = ntrait * sum(nlev)
neq
# make empty lhs and rhs
lhs = matrix(c(0), nrow = neq, ncol = neq)
lhs
rhs = matrix(c(0), nrow = neq)
rhs
# addresses
address = matrix(c(0), nrow=neff, ncol=ntrait)
address
# weight of covariable
weight_cov = matrix(c(0), nrow=neff, ncol=ntrait)
weight_cov
# set working directory setwd(choose.dir())
setwd("D:/temp/00_prediction")
# print working directory
getwd()
# read data : t1e1, t1e2, t2e1, t2e2, y1, y2
data = read.table("06_mt_lsq_data.txt", header = FALSE, sep = "", col.names = c("t1e1", "t1e2", "t2e1", "t2e2", "y1", "y2"), stringsAsFactors = FALSE)
data
# number of data
ndata = nrow(data)
ndata
# read pedigree : animal, sire, dam
pedi = read.table("04_stam_pedi.txt", header = FALSE, sep = "", col.names = c("animal",
"sire", "dam"), stringsAsFactors = FALSE)
pedi
# number of pedigree
nped = nrow(pedi)
nped
# inverse matrix of NRM(Nemerator Relationship Matrix)
a = matrix(c(0), nrow = nped, ncol = nped)
a
# loop data
for (di in 1:ndata) {
# di = 5
indata = as.numeric(data[di, 1:(ntrait * neff)])
y = as.numeric(data[di, (ntrait * neff + 1):(ntrait * (neff + 1))])
# address and weight_cov
find_addresses()
# inverse of r
find_rinv()
for (i in 1:neff) {
for (j in 1:neff) {
for (k in 1:ntrait) {
for (l in 1:ntrait) {
lhs[address[i, k], address[j, l]] = lhs[address[i, k], address[j, l]] + weight_cov[i, k] * weight_cov[j, l] * rinv[k, l]
}
}
}
for (k in 1:ntrait) {
for (l in 1:ntrait) {
rhs[address[i, k]] = rhs[address[i, k]] + weight_cov[i, k] * rinv[k, l] * y[l]
}
}
}
}
# print
lhs
rhs
# random effects' contributions
for(i in 1 : neff){
if(randomtype[i] == 'g_fixed'){
next
} else if(randomtype[i] == 'g_diag'){
add_g_diag(i)
} else if(randomtype[i] == 'g_A' || randomtype[i] == 'g_As'){
add_g_add(randomtype[i], i)
} else {
warning("unimplemented random type")
}
}
# print lhs
lhs
# generalised inverse of lhs
gi_lhs = ginv(lhs)
# print
gi_lhs
# solution
sol = gi_lhs %*% rhs
# print
sol
실행 결과
> # MASS packages
> if (!("MASS" %in% installed.packages())) {
+ install.packages("MASS", dependencies = TRUE)
+ }
> library(MASS)
>
> # function : find address
> find_addresses <- function() {
+
+ for (i in 1:neff) {
+ for (j in 1:ntrait) {
+
+ # i=2 j=1
+
+ if (datum(i, j) == miss) {
+ address[i, j] <<- 0
+ weight_cov[i, j] <<- 0
+ next
+ }
+
+ baseaddr = ifelse(i == 1, j, sum(nlev[1:i - 1]) * ntrait + j)
+
+ if (effecttype[i] == "effcross") {
+ weight_cov[i, j] <<- 1
+ address[i, j] <<- baseaddr + (datum(i, j) - 1) * ntrait
+ } else if (effecttype[i] == "effcov") {
+ weight_cov[i, j] <<- datum(i, j)
+ if (nestedcov[i] == 0) {
+ address[i, j] <<- baseaddr
+ } else if (nestedcov[i] > 0 && nestedcov[i] < neff) {
+ address[i, j] <<- baseaddr + (datum(nestedcov[i], j) - 1) * ntrait
+ } else {
+ warning("wrong description of nested covariable")
+ }
+ } else {
+ warning("unimplemented effect")
+ }
+ }
+ }
+
+ }
>
> # return data in indata
> datum <- function(ef, tr) {
+ return(indata[ef + (tr - 1) * neff])
+ }
>
> # inverse of residual variance and covariance matrix
> find_rinv <- function() {
+ tempr = r
+ for (i in 1:ntrait) {
+ if (y[i] == miss) {
+ tempr[i, ] = 0
+ tempr[, i] = 0
+ }
+ }
+ rinv <<- ginv(tempr)
+ }
>
> # address for given level l of effect e and trait t
> address1 <- function(e, l, t){
+ return(ifelse(e == 1, (l - 1) * ntrait + t, sum(nlev[1 : e - 1]) * ntrait + (l - 1) * ntrait + t))
+ }
>
> # setup g
> setup_g <- function(){
+ for(i in 1 : neff){
+ if (randomnumb[i] != 0){
+ g[,,i] <<- ginv(g0[,,i])
+ }
+ }
+ }
>
> # add a genetic variance to the diagonal of left hand side
> add_g_diag <- function(eff){
+ for(i in 1 : nlev[eff]){
+ for(j in 0 : (randomnumb[eff] - 1)){
+ for(k in 0 : (randomnumb[eff] - 1)){
+ for(t1 in 1 : ntrait){
+ for(t2 in 1 : ntrait){
+ m = address1(eff + j, i, t1)
+ l = address1(eff + k, i, t2)
+ #print(paste("i=",i,"j=",j,"k=",k,"t1=",t1,"t2=",t2,"m=",m,"l=",l))
+ lhs[m, l] <<- lhs[m, l] + g[t1 + j * ntrait, t2 + k * ntrait, eff]
+ }
+ }
+ }
+ }
+ }
+ }
>
> # add genetic variacne matrix * NRN to left hand side
> add_g_add <- function(type, eff){
+
+ if(type == 'g_A'){
+ w = c(1, -0.5, -0.5)
+ res = c(2, 4/3, 1, 0)
+ } else if(type == 'g_As'){
+ w = c(1, -0.5, -0.25)
+ res = c(16/11, 4/3, 16/15, 1)
+ }
+
+ for (di in 1:nped) {
+
+ p = as.numeric(pedi[di,])
+
+ ri = 1
+ if(p[2] == 0){ ri = ri + 1 }
+ if(p[3] == 0){ ri = ri + 1 }
+
+ mendel = res[ri]
+
+ for(i in 0 : (randomnumb[eff] - 1)){
+ for(j in 0 : (randomnumb[eff] - 1)){
+ for(t1 in 1 : ntrait){
+ for(t2 in 1 : ntrait){
+ for (k in 1 : 3) {
+ for (l in 1 : 3) {
+ if (p[k] != 0 && p[l] != 0) {
+ m = address1(eff + i, p[k], t1)
+ n = address1(eff + j, p[l], t2)
+ a[p[k], p[l]] <<- a[p[k], p[l]] + w[k] * w[l] * mendel
+ lhs[m, n] <<- lhs[m, n] + g[t1 + i * ntrait, t2 + j * ntrait, eff] * w[k] * w[l] * mendel
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ }
+
+ }
>
>
> # parameters
>
> # number of traits
> ntrait = 2
> ntrait
[1] 2
>
> # number of effects
> neff = 2
> neff
[1] 2
>
> # type of effects
> effecttype = c("effcross", "effcross")
> effecttype
[1] "effcross" "effcross"
>
> # whether covariable effect is nested or not
> nestedcov = c(0, 0)
> nestedcov
[1] 0 0
>
> # number of levels
> nlev = c(2, 5)
> nlev
[1] 2 5
>
> # missing values
> miss = 0
> miss
[1] 0
>
> # type of random effect : g_fixed, g_diag, g_A, g_As
> randomtype = c("g_fixed", "g_A")
> randomtype
[1] "g_fixed" "g_A"
>
> # number of correlated effects per effect
> randomnumb = c(0, 1)
> randomnumb
[1] 0 1
>
> # genetic varinace and covariance matrix
> g0 = array(c(0), dim =c(ntrait, ntrait, neff))
> g0[,,2] = matrix(c(2, -1, -1, 1), nrow = ntrait)
> g0
, , 1
[,1] [,2]
[1,] 0 0
[2,] 0 0
, , 2
[,1] [,2]
[1,] 2 -1
[2,] -1 1
>
> # matrix for inverse of genetic variance matrix
> g = array(c(0), dim =c(ntrait, ntrait, neff))
> g
, , 1
[,1] [,2]
[1,] 0 0
[2,] 0 0
, , 2
[,1] [,2]
[1,] 0 0
[2,] 0 0
>
> setup_g()
>
> # residual variace
> r = matrix(c(1, 0, 0, 1), nrow = 2)
> r
[,1] [,2]
[1,] 1 0
[2,] 0 1
>
> # empty inverse of r
> rinv = matrix(c(0), nrow = ntrait, ncol = ntrait)
> rinv
[,1] [,2]
[1,] 0 0
[2,] 0 0
>
> # number of equations
> neq = ntrait * sum(nlev)
> neq
[1] 14
>
> # make empty lhs and rhs
> lhs = matrix(c(0), nrow = neq, ncol = neq)
> lhs
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 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
[3,] 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
[5,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[6,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[7,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[8,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[9,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[10,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[11,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[12,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[13,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[14,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
>
> rhs = matrix(c(0), nrow = neq)
> rhs
[,1]
[1,] 0
[2,] 0
[3,] 0
[4,] 0
[5,] 0
[6,] 0
[7,] 0
[8,] 0
[9,] 0
[10,] 0
[11,] 0
[12,] 0
[13,] 0
[14,] 0
>
> # addresses
> address = matrix(c(0), nrow=neff, ncol=ntrait)
> address
[,1] [,2]
[1,] 0 0
[2,] 0 0
>
> # weight of covariable
> weight_cov = matrix(c(0), nrow=neff, ncol=ntrait)
> weight_cov
[,1] [,2]
[1,] 0 0
[2,] 0 0
>
> # set working directory setwd(choose.dir())
> setwd("D:/temp/00_prediction")
>
> # print working directory
> getwd()
[1] "D:/temp/00_prediction"
>
> # read data : t1e1, t1e2, t2e1, t2e2, y1, y2
> data = read.table("06_mt_lsq_data.txt", header = FALSE, sep = "", col.names = c("t1e1", "t1e2", "t2e1", "t2e2", "y1", "y2"), stringsAsFactors = FALSE)
> data
t1e1 t1e2 t2e1 t2e2 y1 y2
1 2 3 2 3 9 1
2 2 2 2 2 6 5
3 2 1 2 1 7 3
4 1 2 1 2 8 4
5 1 1 1 1 5 2
>
> # number of data
> ndata = nrow(data)
> ndata
[1] 5
>
> # read pedigree : animal, sire, dam
> pedi = read.table("04_stam_pedi.txt", header = FALSE, sep = "", col.names = c("animal",
+ "sire", "dam"), stringsAsFactors = FALSE)
> pedi
animal sire dam
1 5 0 0
2 1 3 4
3 2 3 4
4 3 5 0
5 4 0 0
>
> # number of pedigree
> nped = nrow(pedi)
> nped
[1] 5
>
> # inverse matrix of NRM(Nemerator Relationship Matrix)
> a = matrix(c(0), nrow = nped, ncol = nped)
> a
[,1] [,2] [,3] [,4] [,5]
[1,] 0 0 0 0 0
[2,] 0 0 0 0 0
[3,] 0 0 0 0 0
[4,] 0 0 0 0 0
[5,] 0 0 0 0 0
>
> # loop data
> for (di in 1:ndata) {
+ # di = 5
+ indata = as.numeric(data[di, 1:(ntrait * neff)])
+ y = as.numeric(data[di, (ntrait * neff + 1):(ntrait * (neff + 1))])
+
+ # address and weight_cov
+ find_addresses()
+
+ # inverse of r
+ find_rinv()
+
+ for (i in 1:neff) {
+ for (j in 1:neff) {
+ for (k in 1:ntrait) {
+ for (l in 1:ntrait) {
+ lhs[address[i, k], address[j, l]] = lhs[address[i, k], address[j, l]] + weight_cov[i, k] * weight_cov[j, l] * rinv[k, l]
+ }
+ }
+ }
+ for (k in 1:ntrait) {
+ for (l in 1:ntrait) {
+ rhs[address[i, k]] = rhs[address[i, k]] + weight_cov[i, k] * rinv[k, l] * y[l]
+ }
+ }
+ }
+ }
>
> # print
> lhs
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 2 0 0 0 1 0 1 0 0 0 0 0 0 0
[2,] 0 2 0 0 0 1 0 1 0 0 0 0 0 0
[3,] 0 0 3 0 1 0 1 0 1 0 0 0 0 0
[4,] 0 0 0 3 0 1 0 1 0 1 0 0 0 0
[5,] 1 0 1 0 2 0 0 0 0 0 0 0 0 0
[6,] 0 1 0 1 0 2 0 0 0 0 0 0 0 0
[7,] 1 0 1 0 0 0 2 0 0 0 0 0 0 0
[8,] 0 1 0 1 0 0 0 2 0 0 0 0 0 0
[9,] 0 0 1 0 0 0 0 0 1 0 0 0 0 0
[10,] 0 0 0 1 0 0 0 0 0 1 0 0 0 0
[11,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[12,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[13,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[14,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0
> rhs
[,1]
[1,] 13
[2,] 6
[3,] 22
[4,] 9
[5,] 12
[6,] 5
[7,] 14
[8,] 9
[9,] 9
[10,] 1
[11,] 0
[12,] 0
[13,] 0
[14,] 0
>
> # random effects' contributions
> for(i in 1 : neff){
+ if(randomtype[i] == 'g_fixed'){
+ next
+ } else if(randomtype[i] == 'g_diag'){
+ add_g_diag(i)
+ } else if(randomtype[i] == 'g_A' || randomtype[i] == 'g_As'){
+ add_g_add(randomtype[i], i)
+ } else {
+ warning("unimplemented random type")
+ }
+ }
>
> # print lhs
> lhs
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 2 0 0 0 1 0 1 0 0.0000000 0.0000000 0 0 0.0000000 0.0000000
[2,] 0 2 0 0 0 1 0 1 0.0000000 0.0000000 0 0 0.0000000 0.0000000
[3,] 0 0 3 0 1 0 1 0 1.0000000 0.0000000 0 0 0.0000000 0.0000000
[4,] 0 0 0 3 0 1 0 1 0.0000000 1.0000000 0 0 0.0000000 0.0000000
[5,] 1 0 1 0 4 2 0 0 -1.0000000 -1.0000000 -1 -1 0.0000000 0.0000000
[6,] 0 1 0 1 2 6 0 0 -1.0000000 -2.0000000 -1 -2 0.0000000 0.0000000
[7,] 1 0 1 0 0 0 4 2 -1.0000000 -1.0000000 -1 -1 0.0000000 0.0000000
[8,] 0 1 0 1 0 0 2 6 -1.0000000 -2.0000000 -1 -2 0.0000000 0.0000000
[9,] 0 0 1 0 -1 -1 -1 -1 3.3333333 2.3333333 1 1 -0.6666667 -0.6666667
[10,] 0 0 0 1 -1 -2 -1 -2 2.3333333 5.6666667 1 2 -0.6666667 -1.3333333
[11,] 0 0 0 0 -1 -1 -1 -1 1.0000000 1.0000000 2 2 0.0000000 0.0000000
[12,] 0 0 0 0 -1 -2 -1 -2 1.0000000 2.0000000 2 4 0.0000000 0.0000000
[13,] 0 0 0 0 0 0 0 0 -0.6666667 -0.6666667 0 0 1.3333333 1.3333333
[14,] 0 0 0 0 0 0 0 0 -0.6666667 -1.3333333 0 0 1.3333333 2.6666667
>
> # generalised inverse of lhs
> gi_lhs = ginv(lhs)
>
> # print
> gi_lhs
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
[1,] 1.9090909 -0.6969697 1.3333333 -0.6666667 -1.4090909 0.6969697 -1.4090909 0.6969697 -1.1818182 0.6060606 -0.8181818
[2,] -0.6969697 1.2121212 -0.6666667 0.6666667 0.6969697 -0.7121212 0.6969697 -0.7121212 0.6060606 -0.5757576 0.3939394
[3,] 1.3333333 -0.6666667 1.6666667 -0.6666667 -1.3333333 0.6666667 -1.3333333 0.6666667 -1.3333333 0.6666667 -0.6666667
[4,] -0.6666667 0.6666667 -0.6666667 1.0000000 0.6666667 -0.6666667 0.6666667 -0.6666667 0.6666667 -0.6666667 0.3333333
[5,] -1.4090909 0.6969697 -1.3333333 0.6666667 1.5590909 -0.7469697 1.2590909 -0.6469697 1.1818182 -0.6060606 0.8181818
[6,] 0.6969697 -0.7121212 0.6666667 -0.6666667 -0.7469697 0.8121212 -0.6469697 0.6121212 -0.6060606 0.5757576 -0.3939394
[7,] -1.4090909 0.6969697 -1.3333333 0.6666667 1.2590909 -0.6469697 1.5590909 -0.7469697 1.1818182 -0.6060606 0.8181818
[8,] 0.6969697 -0.7121212 0.6666667 -0.6666667 -0.6469697 0.6121212 -0.7469697 0.8121212 -0.6060606 0.5757576 -0.3939394
[9,] -1.1818182 0.6060606 -1.3333333 0.6666667 1.1818182 -0.6060606 1.1818182 -0.6060606 1.6363636 -0.7878788 0.3636364
[10,] 0.6060606 -0.5757576 0.6666667 -0.6666667 -0.6060606 0.5757576 -0.6060606 0.5757576 -0.7878788 0.8484848 -0.2121212
[11,] -0.8181818 0.3939394 -0.6666667 0.3333333 0.8181818 -0.3939394 0.8181818 -0.3939394 0.3636364 -0.2121212 1.6363636
[12,] 0.3939394 -0.4242424 0.3333333 -0.3333333 -0.3939394 0.4242424 -0.3939394 0.4242424 -0.2121212 0.1515152 -0.7878788
[13,] -0.5909091 0.3030303 -0.6666667 0.3333333 0.5909091 -0.3030303 0.5909091 -0.3030303 0.8181818 -0.3939394 0.1818182
[14,] 0.3030303 -0.2878788 0.3333333 -0.3333333 -0.3030303 0.2878788 -0.3030303 0.2878788 -0.3939394 0.4242424 -0.1060606
[,12] [,13] [,14]
[1,] 0.39393939 -0.5909091 0.30303030
[2,] -0.42424242 0.3030303 -0.28787879
[3,] 0.33333333 -0.6666667 0.33333333
[4,] -0.33333333 0.3333333 -0.33333333
[5,] -0.39393939 0.5909091 -0.30303030
[6,] 0.42424242 -0.3030303 0.28787879
[7,] -0.39393939 0.5909091 -0.30303030
[8,] 0.42424242 -0.3030303 0.28787879
[9,] -0.21212121 0.8181818 -0.39393939
[10,] 0.15151515 -0.3939394 0.42424242
[11,] -0.78787879 0.1818182 -0.10606061
[12,] 0.84848485 -0.1060606 0.07575758
[13,] -0.10606061 1.9090909 -0.94696970
[14,] 0.07575758 -0.9469697 0.96212121
>
> # solution
> sol = gi_lhs %*% rhs
>
> # print
> sol
[,1]
[1,] 7.0606061
[2,] 2.5757576
[3,] 7.3333333
[4,] 3.0000000
[5,] -0.6606061
[6,] 0.1242424
[7,] -0.4606061
[8,] 0.7242424
[9,] 1.1212121
[10,] -0.8484848
[11,] -1.1212121
[12,] 0.8484848
[13,] 0.5606061
[14,] -0.4242424
-
기본적인 가축 육종 모형의 육종가를 구하는 프로그램
복잡한 것은 모르겠고, 다형질 임의회귀 모형(multiple trait random regression model)의 육종가까지는 구할 수 있다. 임의회귀 모형에서는 개체별로 회귀식을 추정하므로(개체별로 nested된 회귀식) MME에는 들어가지 않지만 어느 개체에 nested 할지를 알려 주어야 하므로 그에 따른 처리가 필요하다.
프로그램
# Best Linear Unbiase Prediction of Breeding Values
# Date : 2020.08.06.
# Reference : Computational techneques in animal breeding, Ignacy Misztal
# MASS packages
if (!("MASS" %in% installed.packages())) {
install.packages("MASS", dependencies = TRUE)
}
library(MASS)
# function : find address
find_addresses <- function() {
for (i in 1:neff) {
for (j in 1:ntrait) {
if (datum(i, j) == miss) {
address[i, j] <<- 0
weight_cov[i, j] <<- 0
next
}
baseaddr = ifelse(i == 1, j, sum(nlev[1:i - 1]) * ntrait + j)
if (effecttype[i] == "effcross") {
weight_cov[i, j] <<- 1
address[i, j] <<- baseaddr + (datum(i, j) - 1) * ntrait
} else if (effecttype[i] == "effcov") {
weight_cov[i, j] <<- datum(i, j)
if (nestedcov[i] == 0) {
address[i, j] <<- baseaddr
} else if (nestedcov[i] > 0 && nestedcov[i] < neff) {
address[i, j] <<- baseaddr + (datum(nestedcov[i], j) - 1) * ntrait
} else {
warning("wrong description of nested covariable")
}
} else if(effecttype[i] == "effnone") {
} else {
warning("unimplemented effect")
}
}
}
}
# return data in indata
datum <- function(ef, tr) {
return(indata[ef + (tr - 1) * neff])
}
# inverse of residual variance and covariance matrix
find_rinv <- function() {
tempr = r
for (i in 1:ntrait) {
if (y[i] == miss) {
tempr[i, ] = 0
tempr[, i] = 0
}
}
rinv <<- ginv(tempr)
}
# address for given level l of effect e and trait t
address1 <- function(e, l, t){
return(ifelse(e == 1, (l - 1) * ntrait + t, sum(nlev[1 : e - 1]) * ntrait + (l - 1) * ntrait + t))
}
# setup g
setup_g <- function(){
for(i in 1 : neff){
if (randomnumb[i] != 0){
g[,,i] <<- ginv(g0[,,i])
}
}
}
# add a genetic variance to the diagonal of left hand side
add_g_diag <- function(eff){
for(i in 1 : nlev[eff]){
for(j in 0 : (randomnumb[eff] - 1)){
for(k in 0 : (randomnumb[eff] - 1)){
for(t1 in 1 : ntrait){
for(t2 in 1 : ntrait){
m = address1(eff + j, i, t1)
l = address1(eff + k, i, t2)
#print(paste("i=",i,"j=",j,"k=",k,"t1=",t1,"t2=",t2,"m=",m,"l=",l))
lhs[m, l] <<- lhs[m, l] + g[t1 + j * ntrait, t2 + k * ntrait, eff]
}
}
}
}
}
}
# add genetic variacne matrix * NRN to left hand side
add_g_add <- function(type, eff){
if(type == 'g_A'){
w = c(1, -0.5, -0.5)
res = c(2, 4/3, 1, 0)
} else if(type == 'g_As'){
w = c(1, -0.5, -0.25)
res = c(16/11, 4/3, 16/15, 1)
}
for (di in 1:nped) {
p = as.numeric(pedi[di,])
ri = 1
if(p[2] == 0){ ri = ri + 1 }
if(p[3] == 0){ ri = ri + 1 }
mendel = res[ri]
for(i in 0 : (randomnumb[eff] - 1)){
for(j in 0 : (randomnumb[eff] - 1)){
for(t1 in 1 : ntrait){
for(t2 in 1 : ntrait){
for (k in 1 : 3) {
for (l in 1 : 3) {
if (p[k] != 0 && p[l] != 0) {
m = address1(eff + i, p[k], t1)
n = address1(eff + j, p[l], t2)
a[p[k], p[l]] <<- a[p[k], p[l]] + w[k] * w[l] * mendel
lhs[m, n] <<- lhs[m, n] + g[t1 + i * ntrait, t2 + j * ntrait, eff] * w[k] * w[l] * mendel
}
}
}
}
}
}
}
}
}
# parameters
# number of traits
ntrait = 2
ntrait
# number of effects
neff = 2
neff
# type of effects : cross-classified or covariable
effecttype = c("effcross", "effcross")
effecttype
# 0 means cross-classified or not nested, 1 means # of effect within which covariable is nested
nestedcov = c(0, 0)
nestedcov
# level of effect, level of covariable depends on if it is nested or not
nlev = c(2, 5)
nlev
# missing values trait/effect
miss = 0
miss
# type of random effect : g_fixed, g_diag, g_A, g_As
randomtype = c("g_fixed", "g_A")
randomtype
# number of correlated effects per effect
randomnumb = c(0, 1)
randomnumb
# genetic varinace and covariance matrix for each trait
vardim = ntrait * max(randomnumb)
g0 = array(c(0), dim =c(vardim, vardim, neff))
# genetic variance matrix for effect 2
g0[,,2] = matrix(c(2, -1, -1, 1), nrow = vardim)
g0
# matrix for inverse of genetic variance matrix
g = array(c(0), dim =c(vardim, vardim, neff))
setup_g()
g
# residual variace
r = matrix(c(1, 0, 0, 1), nrow = ntrait)
r
# empty inverse of r
rinv = matrix(c(0), nrow = ntrait, ncol = ntrait)
rinv
# number of equations
neq = ntrait * sum(nlev)
neq
# make empty lhs and rhs
lhs = matrix(c(0), nrow = neq, ncol = neq)
lhs
rhs = matrix(c(0), nrow = neq)
rhs
# addresses
address = matrix(c(0), nrow=neff, ncol=ntrait)
address
# weight of covariable
weight_cov = matrix(c(0), nrow=neff, ncol=ntrait)
weight_cov
# set working directory setwd(choose.dir())
setwd("D:/temp/00_prediction")
# print working directory
getwd()
# read data : t1e1, t1e2, t2e1, t2e2, y1, y2
data = read.table("06_mt_lsq_data.txt", header = FALSE, sep = "", col.names = c("t1e1", "t1e2", "t2e1", "t2e2", "y1", "y2"), stringsAsFactors = FALSE)
data
# number of data
ndata = nrow(data)
ndata
# read pedigree : animal, sire, dam
pedi = read.table("04_stam_pedi.txt", header = FALSE, sep = "", col.names = c("animal",
"sire", "dam"), stringsAsFactors = FALSE)
pedi
# number of pedigree
nped = nrow(pedi)
nped
# inverse matrix of NRM(Nemerator Relationship Matrix)
a = matrix(c(0), nrow = nped, ncol = nped)
a
# loop data
for (di in 1:ndata) {
indata = as.numeric(data[di, 1:(ntrait * neff)])
y = as.numeric(data[di, (ntrait * neff + 1):(ntrait * (neff + 1))])
# address and weight_cov
find_addresses()
# inverse of r
find_rinv()
for (i in 1:neff) {
for (j in 1:neff) {
for (k in 1:ntrait) {
for (l in 1:ntrait) {
lhs[address[i, k], address[j, l]] = lhs[address[i, k], address[j, l]] + weight_cov[i, k] * weight_cov[j, l] * rinv[k, l]
}
}
}
for (k in 1:ntrait) {
for (l in 1:ntrait) {
rhs[address[i, k]] = rhs[address[i, k]] + weight_cov[i, k] * rinv[k, l] * y[l]
}
}
}
}
# print
lhs
rhs
# random effects' contributions
for(i in 1 : neff){
if(randomtype[i] == 'g_fixed'){
next
} else if(randomtype[i] == 'g_diag'){
add_g_diag(i)
} else if(randomtype[i] == 'g_A' || randomtype[i] == 'g_As'){
add_g_add(randomtype[i], i)
} else {
warning("unimplemented random type")
}
}
# print lhs
lhs
# generalised inverse of lhs
gi_lhs = ginv(lhs)
# print
gi_lhs
# solution
sol = gi_lhs %*% rhs
# print
sol
# write a solution by effect, level and trait
# copy the continuous variances to the original effects
for(i in 1 : neff){
if (randomnumb[i] > 1){
for(j in 2 : randomnumb[i]){
first = ntrait * (j - 1) + 1
last = ntrait * j
print(paste(i, j, first, last, (neff + j - 1)))
g0[1 : ntrait, 1 : ntrait, (i + j - 1)] = g0[first : last, first : last, i]
}
}
}
line = c("effect level trait solution sep rel")
write(line, file = "solutions")
for(e in 1 : neff){
# when the effect is used
if(nlev[e] > 0){
for(l in 1 : nlev[e]){
for(t in 1 : ntrait){
j = address1(e, l, t)
sep2 = gi_lhs[j, j]
# when the effect is random
rel = ifelse(randomtype[e] != 'g_fixed',1 - sep2 / g0[t, t, e], 0)
line = paste(e, l, t, sol[j], sqrt(sep2), rel)
write(line, file = "solutions", append = TRUE)
}
}
}
}
마지막 solution을 출력하는 부분은 개선이 필요하다.
'Animal Breeding > R for Genetic Evaluation' 카테고리의 다른 글
| 임의 회귀 모형의 해를 이용한 개체 육종가 그래프 및 305일령 육종가 계산 (0) | 2020.11.26 |
|---|---|
| 다형질 모형의 차원 줄이기(Reduce the Dimension of Multivariate Models) (0) | 2020.09.29 |
| R을 이용한 다형질 개체 모형(unequal design matrices)의 육종가 구하기 프로그래밍 (0) | 2020.07.29 |
| R을 이용한 다형질 개체 모형(same model and missing record)의 육종가 구하기 프로그래밍 (0) | 2020.07.26 |
| R을 이용한 다형질 개체 모형(same model and no missing record)의 육종가 구하기 프로그래밍 (0) | 2020.07.03 |