레노보 아이디어패드 스림3-15IIL05

LENOVO IdeaPad Slim3-15 IIL05

 

또 노트북을 구매했다. 옛날에 라디오 들으며 공부했다면 요즘 아이들은 유투브로 음악을 들으면서 공부를 하는가 보다. 그러다 보니 컴퓨터가 내 차례가 되질 않는다. 그래서 질렀다. 눈이 침침하니 크면 좋겠으나 지난번에 산 17인치는 사실 들고 다니는 노트북이 아니다. 그리고 키패드가 있으면 좋겠다. 14인치 이하는 키패드가 없다. 결국 15인치밖에 살게 없다. 역시 Lenovo IdeaPad Slim3가 저렴하다. 지난번에 산 17인치하고 다른 점은 화면이 크고, 배터리 용량이 3 cell(15인치는 2 cell) 이란 거 외에 차이점은 없는 거 같다. 인텔 Core i3이지만 10세대라 부팅속도나 크롬 브라우저 시작 속도가 빠르다. 그냥 간단한 인터넷, 영화보기에 이 정도면 부족함이 없어 보인다.

 

한 달에 노트북을 2대나 샀다. 안사람까지 1인 1노트북이 되었다.

 

 

 

박스 앞 면

 

 

박스 뒷면

 

 

본체와 충전기 포장 상태

 

 

본체와 충전기의 포장을 벗김. HDD를 달 수 있는 키트와 그냥 별 내용 없는 하지만 제조자 입장에서는 안 넣으면 안 되는 종이들.

 

 

위에 바라본 노트북. 그리 싼 티는 나지않는다. 플라스틱처럼 보이지는 않는다. 뭔지는 정확히 잘 모르겠다.

 

 

왼쪽. 전원, HDMI, 세 개의 USB

 

 

노트북 오른쪽. 이어폰 구멍. 마이크. SD 카드 리더.

 

 

뚜껑을 연 모습.

 

 

보호 종이 제거.

 

 

뒷 면은 플라스틱이다.

 

 

램과 SSD를 추가하기 위하여 뒷판을 열었다.

 

 

램을 추가하고, 가지고 있던 512G SSD 설치

 

 

8G 램 추가

 

 

원래 달려 있던 256G NVMe SSD. 웨스턴 디지탈 것인가 보다.

 

 

HDMI, USB 등

 

 

이어폰하고 SD 카드 리더.

 

 

뒷판을 덮었다.

 

 

 

이해가 가지 않는 것은 하드를 달기 어럽게 만들어 놓았다는 점이다. 배터리 고정 나사까지 풀고, HDD 업그레이드 키트의 선을 이용하여 메인보드와 연결하였다. 그냥 특정 부위에 하드를 끼면 되게 만들면 될 것을. 아마도 원가 절감 차원에서 그런 거 같다. 42만원 안 되게 주었으니 그 정도의 불편은 참는 수밖에. 

인텔 Core i3 10th Gen. RAM 12Gb, NVMe 256G SSD, 512SSD. 이 정도면 good.

 

윈도우는 포함안되어 있음. 40만원대 컴퓨터에 너무 많은 걸 바라면 안된다.

 

Lenovo IdeaPad Slim3-17IML 3D 구입

만 8년이 된 ASUS 노트북은 아직도 고장이 안나고 잘 돌아간다. ASUS 노트북은 죽지 않는다. 다만 느려질 뿐이다. 온라인 클래스만 아니어도 안 바꾸는데 온라인 클래스를 느린 컴퓨터로 듣는게 힘들다는 아이의 말에 바꾸지 않을 수가 없었다. 사실은 게임 관련 유투브 볼 때 고해상도로 보고 싶은 것이었는데 ...

느린 컴퓨터가 나쁜 것만은 아닌다. 아이가 게임을 하지 못한다. League of Legends를 하고 싶어도 안 돌아간다. 그래서 새 컴퓨터를 살 때의 기준은 인터넷, 온라인 클래스는 어느 정도 되지만 게임은 안되는 정도의 PC가 필요했다. 즉 현 시점에서 약간 후진 컴퓨터가 필요했다.(사실은 돈이 없어 저렴한 컴퓨터가 필요했다.)

데스크톱은 부피가 커서 이번에도 노트북 중 가장 화면이 큰 17인치 중에서, 그리고 적당한 CPU(Core i3)로 골랐다. 그래도 몇 년 쓸 노트북인데 셀러론이나 펜티엄 골드는 피하고 싶었다. 집에 놓고 항상 전기를 연결할 수 있으므로 배터리는 중요하지 않았다.

그리하여 고른 것이 Lenovo IdeaPad Slim3-17IML 3D 이다.

 

박스

 

 

 

노트북

 

 

 

왼쪽면. 전원. HDMI, USB 3개. 오른쪽 2개 USB가 고속.

 

 

 

오른쪽. 헤드폰, 마이크. SD 카드 리더. 이전 노트북(ASUS 노트북)은 CD-ROM 드라이브가 있었는데 이번 모델은 없다. 당장 CD로 설치하는 프로그램이 있는데 깔 수가 없다. 외장 CD-ROM 드라이브를 빌려서 깔아야 한다.

 

 

 

펼쳤을 때. 카메라 셔터가 있어서 안 쓸 때는 닫을 수 있다. 

 

 

 

뒷면. 뒷면을 열고 8G 램을 설치하고, 1T HDD를 추가 했다. SLIM3가 단점도 있지만 이렇게 램과 HDD를 추가할 수 있는 점은 장점이라고 할 수 있다. 

 

 

 

윈도우즈를 설치하고 부팅한 화면. 부팅은 정말 빠르다. 10초 이내인 것 같다. 그리고 노트북 뚜껑을 열면 바로 켜진다.

 

 

 

인터넷, 온라이 클래스, 문서, 엑셀, PPT 작업, 간단한 영화 보기에 적당한 컴퓨터이다. 이걸로 엄청난 게임을 못할 것이고, 동영상 편집도 못할 것이지만 난 그런 거 할 계획이 없으므로 우리집에 알맞은 컴퓨터로 보인다. 

 

전체적인 기능은 비슷한데 자잘한 기능을 자꾸 넣어서 OS를 포함한 소프트웨어가 자꾸 느려지니 하드웨어를 바꿀 수 밖에 없다. 

이건 몇 년 쓸 수 있을까?

다형질 모형을 분석하는 것은 매우 큰 행렬을 다루는 것이다. 예를 1000마리를 유전평가한다고 하였을 때 single trait model의 경우 MME의 원소 수는 1000 * 1000으로 백만 개의 원소이나, 2 trait model의 경우 2000 * 2000으로 원소 수가 4백만 개로 늘어난다. n traits가 되면 원소 수는 n의 제곱이 된다. 그래서 다형질 모형을 자료 변환을 통하여 단형질 모형으로 분석할 방법을 고안하게 되었다.

 

여기서는 Canonical transformation을 통해서 다형질 모형을 단형질 모형으로 분석하는 방법을 소개한다. canonical transformaion matrix을 만들어 원래 자료에서 변환 자료를 만들고, 각각의 변환 자료를 단형질 모형으로 분석한다. 이후 단형질 모형의 육종가를 다시 canonical transformation matrix의 역행렬을 이용하여 변환하면 다형질 모형의 육종가와 동일하게 된다. 

 

단 이 방법은 각 형질의 모형이 동일한 경우(same model) 그리고 결측치가 없는 경우(no missing value)에 사용할 수 있어 실제적으로 사용하기에 상당한 제한이 있다.

 

Linear Models for the Prediction of Animal Breeding Values 3rd Edition, Raphael A. Mrode. 

Example 6.1 p97.

Example 5.1 p72에 나오는 다형질 모형을 canonical transformation matrix 이용하여 단형질 모형으로 분석하고 다시 단형질 분석의 육종가를 변환하는 과정

 

자료 mt01_data.txt

4 1 4.5 6.8
5 2 2.9 5.0
6 2 3.9 6.8
7 1 3.5 6.0
8 1 5.0 7.5

 

혈통 mt01_pedi.txt

1 0 0
2 0 0
3 0 0
4 1 0
5 3 2
6 1 2
7 4 5
8 3 6

 

프로그램

# Reduce the dimension of multivariate models
# canonical transformation 
# n traits model to single trait model

# MASS packages
if (!("MASS" %in% installed.packages())) {
    install.packages("MASS", dependencies = TRUE)
}
library(MASS)

# plyr packages
if (!("plyr" %in% installed.packages())) {
    install.packages("plyr", dependencies = TRUE)
}
library(plyr)

# genetic variance - covariance matrix
g = matrix(c(20, 18, 18, 40), nrow = 2)
g

# residual matrix
r = matrix(c(40, 11, 11, 30), nrow = 2)
r

# eigenvalues and eigenvectors
e1 = eigen(r)
e1

# diagonal matirx of eigenvalues
b = diag(e1$values)
b

# eigenvectors
u = e1$vectors
u

# calculate p
p = u %*% sqrt(ginv(b)) %*% t(u)
p

# pgp
pgp = p %*% g %*% t(p)
pgp

# eigenvalues and eigenvectors
e2 = eigen(pgp)
e2

# eigenvectors
l = e2$vectors
l

# canonical transformation matrix
q = t(l) %*% p
q

# new genetic variance
w = q %*% g %*% t(q)
w

# set working directory 
setwd(".") 

# print working directory 
getwd()

# read data : animal, dam, sex, weaning_weight
data = read.table("mt01_data.txt", header = FALSE, sep = "", col.names = c("animal", "sex", "wwg", "pwg"), stringsAsFactors = FALSE)

# print
data

# number of data
ndata = nrow(data)
ndata

# data2
data2 = cbind(data, wwg2 = rep(0, 5), pwg2 = rep(0, 5))
data2

q %*% as.numeric(data2[1, 3:4])

# data transformation
for (i in 1 : ndata) {
    data2[i, 5:6] = q %*% as.numeric(data2[i, 3:4])
}

data2

# wwg2 data
wwg2 = data2[, c('animal','sex','wwg2')]
wwg2

#pwg2
pwg2 = data2[, c('animal','sex','pwg2')]
pwg2

# output wwg2
write.table(wwg2, 'wwg2_data.txt', sep = " ", row.names = FALSE, col.names = FALSE, quote = FALSE)

# output pwg2
write.table(pwg2, 'pwg2_data.txt', sep = " ", row.names = FALSE, col.names = FALSE, quote = FALSE)

#########################################################
# breeding values estimation using renumf90 and blupf90 #
#########################################################

# read solutions
sol_wwg2 = read.table("solutions_wwg2", skip = 1, sep = "", col.names = c("trait", "effect", "level", "solution", "se"), stringsAsFactors = FALSE)
sol_wwg2

# read solutions
sol_pwg2 = read.table("solutions_pwg2", skip = 1, sep = "", col.names = c("trait", "effect", "level", "solution", "se"), stringsAsFactors = FALSE)
sol_pwg2

# breeding values of wwg2
bv_wwg2 = sol_wwg2[sol_wwg2$effect == 2, c(3, 4)]
bv_wwg2

# breeding values of wwg2
bv_pwg2 = sol_pwg2[sol_pwg2$effect == 2, 4]
bv_pwg2

# breeding values
bv_01 = data.frame(animal = bv_wwg2$level, wwg2 = bv_wwg2$solution, pwg2 = bv_pwg2)
bv_01

# read renumbered pedigree 
pedi = read.table("renadd02.ped", sep = "", stringsAsFactors = FALSE) 
pedi

aid_renum = pedi[,c(1,10)]
colnames(aid_renum) = c("animal", "orig_aid")
aid_renum

# join original animal id
bv_02 = join(x = bv_01, y = aid_renum, by = "animal", type = "left")
bv_02

# sort by original id
bv_02 = bv_02[order(bv_02$orig_aid),]
bv_02

# number of animals in pedigree
npedi = nrow(bv_02)
npedi

# add wwg and pwg
bv_03 = cbind(bv_02, wwg = rep(0, npedi), pwg = rep(0, npedi))
bv_03

# inverse of q
iq = ginv(q)
iq

for (i in 1 : npedi) {
    bv_03[i, 5:6] = iq %*% as.numeric(bv_03[i, 2:3])
}
bv_03

 

실행 결과

> # Reduce the dimension of multivariate models
> # canonical transformation 
> # n traits model to single trait model
> 
> # MASS packages
> if (!("MASS" %in% installed.packages())) {
+     install.packages("MASS", dependencies = TRUE)
+ }
> library(MASS)
> 
> # plyr packages
> if (!("plyr" %in% installed.packages())) {
+     install.packages("plyr", dependencies = TRUE)
+ }
> library(plyr)
> 
> # genetic variance - covariance matrix
> g = matrix(c(20, 18, 18, 40), nrow = 2)
> g
     [,1] [,2]
[1,]   20   18
[2,]   18   40
> 
> # residual matrix
> r = matrix(c(40, 11, 11, 30), nrow = 2)
> r
     [,1] [,2]
[1,]   40   11
[2,]   11   30
> 
> # eigenvalues and eigenvectors
> e1 = eigen(r)
> e1
eigen() decomposition
$values
[1] 47.08305 22.91695

$vectors
           [,1]       [,2]
[1,] -0.8407743  0.5413857
[2,] -0.5413857 -0.8407743

> 
> # diagonal matirx of eigenvalues
> b = diag(e1$values)
> b
         [,1]     [,2]
[1,] 47.08305  0.00000
[2,]  0.00000 22.91695
> 
> # eigenvectors
> u = e1$vectors
> u
           [,1]       [,2]
[1,] -0.8407743  0.5413857
[2,] -0.5413857 -0.8407743
> 
> # calculate p
> p = u %*% sqrt(ginv(b)) %*% t(u)
> p
            [,1]        [,2]
[1,]  0.16424710 -0.02874736
[2,] -0.02874736  0.19038107
> 
> # pgp
> pgp = p %*% g %*% t(p)
> pgp
          [,1]      [,2]
[1,] 0.4026185 0.2643755
[2,] 0.2643755 1.2692999
> 
> # eigenvalues and eigenvectors
> e2 = eigen(pgp)
> e2
eigen() decomposition
$values
[1] 1.3435798 0.3283387

$vectors
          [,1]       [,2]
[1,] 0.2704897 -0.9627229
[2,] 0.9627229  0.2704897

> 
> # eigenvectors
> l = e2$vectors
> l
          [,1]       [,2]
[1,] 0.2704897 -0.9627229
[2,] 0.9627229  0.2704897
> 
> # canonical transformation matrix
> q = t(l) %*% p
> q
            [,1]       [,2]
[1,]  0.01675141 0.17550834
[2,] -0.16590031 0.07917187
> 
> # new genetic variance
> w = q %*% g %*% t(q)
> w
              [,1]          [,2]
[1,]  1.343580e+00 -1.110223e-16
[2,] -4.163336e-17  3.283387e-01
> 
> # set working directory 
> setwd(".") 
> 
> # print working directory 
> getwd()
[1] "D:/temp/07_multi2single"
> 
> # read data : animal, dam, sex, weaning_weight
> data = read.table("mt01_data.txt", header = FALSE, sep = "", col.names = c("animal", "sex", "wwg", "pwg"), stringsAsFactors = FALSE)
> 
> # print
> data
  animal sex wwg pwg
1      4   1 4.5 6.8
2      5   2 2.9 5.0
3      6   2 3.9 6.8
4      7   1 3.5 6.0
5      8   1 5.0 7.5
> 
> # number of data
> ndata = nrow(data)
> ndata
[1] 5
> 
> # data2
> data2 = cbind(data, wwg2 = rep(0, 5), pwg2 = rep(0, 5))
> data2
  animal sex wwg pwg wwg2 pwg2
1      4   1 4.5 6.8    0    0
2      5   2 2.9 5.0    0    0
3      6   2 3.9 6.8    0    0
4      7   1 3.5 6.0    0    0
5      8   1 5.0 7.5    0    0
> 
> q %*% as.numeric(data2[1, 3:4])
           [,1]
[1,]  1.2688381
[2,] -0.2081827
> 
> # data transformation
> for (i in 1 : ndata) {
+     data2[i, 5:6] = q %*% as.numeric(data2[i, 3:4])
+ }
> 
> data2
  animal sex wwg pwg      wwg2        pwg2
1      4   1 4.5 6.8 1.2688381 -0.20818267
2      5   2 2.9 5.0 0.9261208 -0.08525154
3      6   2 3.9 6.8 1.2587872 -0.10864249
4      7   1 3.5 6.0 1.1116800 -0.10561986
5      8   1 5.0 7.5 1.4000696 -0.23571251
> 
> # wwg2 data
> wwg2 = data2[, c('animal','sex','wwg2')]
> wwg2
  animal sex      wwg2
1      4   1 1.2688381
2      5   2 0.9261208
3      6   2 1.2587872
4      7   1 1.1116800
5      8   1 1.4000696
> 
> #pwg2
> pwg2 = data2[, c('animal','sex','pwg2')]
> pwg2
  animal sex        pwg2
1      4   1 -0.20818267
2      5   2 -0.08525154
3      6   2 -0.10864249
4      7   1 -0.10561986
5      8   1 -0.23571251
> 
> # output wwg2
> write.table(wwg2, 'wwg2_data.txt', sep = " ", row.names = FALSE, col.names = FALSE, quote = FALSE)
> 
> # output pwg2
> write.table(pwg2, 'pwg2_data.txt', sep = " ", row.names = FALSE, col.names = FALSE, quote = FALSE)
> 
> #########################################################
> # breeding values estimation using renumf90 and blupf90 #
> #########################################################
> 
> # read solutions
> sol_wwg2 = read.table("solutions_wwg2", skip = 1, sep = "", col.names = c("trait", "effect", "level", "solution", "se"), stringsAsFactors = FALSE)
> sol_wwg2
   trait effect level    solution        se
1      1      1     1  1.26648926 1.0167961
2      1      1     2  1.08895002 1.1351950
3      1      2     1 -0.08935869 1.0261539
4      1      2     2 -0.00239531 0.9989446
5      1      2     3  0.09539732 1.0284425
6      1      2     4  0.07287388 0.9893837
7      1      2     5 -0.08838927 1.0012954
8      1      2     6  0.05159845 1.0980228
9      1      2     7 -0.00159342 1.1410714
10     1      2     8 -0.03120867 1.0653321
> 
> # read solutions
> sol_pwg2 = read.table("solutions_pwg2", skip = 1, sep = "", col.names = c("trait", "effect", "level", "solution", "se"), stringsAsFactors = FALSE)
> sol_pwg2
   trait effect level    solution        se
1      1      1     1 -0.18510837 0.7112092
2      1      1     2 -0.09805258 0.8369194
3      1      2     1  0.01452058 0.5482942
4      1      2     2  0.00069546 0.5417348
5      1      2     3 -0.00480559 0.5484205
6      1      2     4 -0.00940595 0.5391831
7      1      2     5  0.00701672 0.5418626
8      1      2     6 -0.00290040 0.5606925
9      1      2     7  0.00195090 0.5697311
10     1      2     8 -0.00048094 0.5543938
> 
> # breeding values of wwg2
> bv_wwg2 = sol_wwg2[sol_wwg2$effect == 2, c(3, 4)]
> bv_wwg2
   level    solution
3      1 -0.08935869
4      2 -0.00239531
5      3  0.09539732
6      4  0.07287388
7      5 -0.08838927
8      6  0.05159845
9      7 -0.00159342
10     8 -0.03120867
> 
> # breeding values of wwg2
> bv_pwg2 = sol_pwg2[sol_pwg2$effect == 2, 4]
> bv_pwg2
[1]  0.01452058  0.00069546 -0.00480559 -0.00940595  0.00701672 -0.00290040  0.00195090 -0.00048094
> 
> # breeding values
> bv_01 = data.frame(animal = bv_wwg2$level, wwg2 = bv_wwg2$solution, pwg2 = bv_pwg2)
> bv_01
  animal        wwg2        pwg2
1      1 -0.08935869  0.01452058
2      2 -0.00239531  0.00069546
3      3  0.09539732 -0.00480559
4      4  0.07287388 -0.00940595
5      5 -0.08838927  0.00701672
6      6  0.05159845 -0.00290040
7      7 -0.00159342  0.00195090
8      8 -0.03120867 -0.00048094
> 
> # read renumbered pedigree 
> pedi = read.table("renadd02.ped", sep = "", stringsAsFactors = FALSE) 
> pedi
  V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
1  1  2  5  1  0  2  1  0  0   7
2  7  0  0  3  0  0  0  0  2   2
3  2  6  0  2  0  1  1  1  0   4
4  3  6  7  1  0  2  1  0  1   6
5  6  0  0  3  0  0  0  2  0   1
6  4  8  3  1  0  2  1  0  0   8
7  8  0  0  3  0  0  0  2  0   3
8  5  8  7  1  0  2  1  0  1   5
> 
> aid_renum = pedi[,c(1,10)]
> colnames(aid_renum) = c("animal", "orig_aid")
> aid_renum
  animal orig_aid
1      1        7
2      7        2
3      2        4
4      3        6
5      6        1
6      4        8
7      8        3
8      5        5
> 
> # join original animal id
> bv_02 = join(x = bv_01, y = aid_renum, by = "animal", type = "left")
> bv_02
  animal        wwg2        pwg2 orig_aid
1      1 -0.08935869  0.01452058        7
2      2 -0.00239531  0.00069546        4
3      3  0.09539732 -0.00480559        6
4      4  0.07287388 -0.00940595        8
5      5 -0.08838927  0.00701672        5
6      6  0.05159845 -0.00290040        1
7      7 -0.00159342  0.00195090        2
8      8 -0.03120867 -0.00048094        3
> 
> # sort by original id
> bv_02 = bv_02[order(bv_02$orig_aid),]
> bv_02
  animal        wwg2        pwg2 orig_aid
6      6  0.05159845 -0.00290040        1
7      7 -0.00159342  0.00195090        2
8      8 -0.03120867 -0.00048094        3
2      2 -0.00239531  0.00069546        4
5      5 -0.08838927  0.00701672        5
3      3  0.09539732 -0.00480559        6
1      1 -0.08935869  0.01452058        7
4      4  0.07287388 -0.00940595        8
> 
> # number of animals in pedigree
> npedi = nrow(bv_02)
> npedi
[1] 8
> 
> # add wwg and pwg
> bv_03 = cbind(bv_02, wwg = rep(0, npedi), pwg = rep(0, npedi))
> bv_03
  animal        wwg2        pwg2 orig_aid wwg pwg
6      6  0.05159845 -0.00290040        1   0   0
7      7 -0.00159342  0.00195090        2   0   0
8      8 -0.03120867 -0.00048094        3   0   0
2      2 -0.00239531  0.00069546        4   0   0
5      5 -0.08838927  0.00701672        5   0   0
3      3  0.09539732 -0.00480559        6   0   0
1      1 -0.08935869  0.01452058        7   0   0
4      4  0.07287388 -0.00940595        8   0   0
> 
> # inverse of q
> iq = ginv(q)
> iq
         [,1]       [,2]
[1,] 2.600648 -5.7651217
[2,] 5.449516  0.5502527
> 
> for (i in 1 : npedi) {
+     bv_03[i, 5:6] = iq %*% as.numeric(bv_03[i, 2:3])
+ }
> bv_03
  animal        wwg2        pwg2 orig_aid         wwg          pwg
6      6  0.05159845 -0.00290040        1  0.15091058  0.279590613
7      7 -0.00159342  0.00195090        2 -0.01539110 -0.007609879
8      8 -0.03120867 -0.00048094        3 -0.07839010 -0.170336777
2      2 -0.00239531  0.00069546        4 -0.01023877 -0.012670601
5      5 -0.08838927  0.00701672        5 -0.27032165 -0.477817750
3      3  0.09539732 -0.00480559        6  0.27579969  0.517224909
1      1 -0.08935869  0.01452058        7 -0.31610344 -0.478971601
4      4  0.07287388 -0.00940595        8  0.24374578  0.391951708
> 

 

중간에 eigenvector가 책의 결과와 다르다. 그러나 eigenvector는 유일한 것이 아니다. 그냥 무시하고 진행하고 결과를 비교해 보니 동일하였다.

지금까지 여러 가지 모형의 육종가를 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을 출력하는 부분은 개선이 필요하다.

 

앞에서는 각 형질의 모형이 같다고 가정하였다. 즉 각 형질에 동일한 고정 효과와 임의 효과(effect)가 들어간다는 뜻이다. 각 형질의 모형이 다를 때, 즉 모형에 포함된 효과(effect)가 다른 경우를 다룬다. 예를 들어 hys1은 첫번째 형질에만 해당되고, hys2는 두번째 형질에만 해당되는 경우이다. 비록 모형이 다르더라도 모형이 같다고 가정해야 traits in animal(effect) 방식으로 프로그램을 할 수 있다. 

예를 들어 자료가 다음과 같다고 가정하자.

4 1 1 201 280
5 1 2 150 200
6 2 1 160 190
7 1 1 180 250
8 2 2 285 300

 

첫째 컬럼은 animal이다. 둘째 컬럼은 hys1인데 첫째 형질의 모형에만 들어간다. 셋째 컬럼은 hys2인데 둘째 형질의 모형에만 들어간다. 넷째는 fat1으로 1산 유지량 형질, 다섯째는 fat2로 2산 유지량이다. 비록 hys1이 첫째 형질에만 hys2가 둘째 형질에만 있지만, 두 고정효과 모두 양 쪽 형질의 모형에 있다고 가정한다. 그래서 same model이라고 하면 다음과 같이 design matrix가 만들어질 것이다.

 

그러나 사실 hys1은 형질 2에 없고, hys2는 형질 1에 없으므로 다음과 같이 만들어진다.

 

위는 효과 내에, 형질 내에, 레벨이 나오는데 이것은 형질 -> 레벨 -> 형질의 순으로 바꾸어 생각한다.

위의 design matrix에서 4번 개체가 LHS에 기여하는 바는 다음과 같다. R의 inverse를 곱하여 얻은 값이다.

 

총 9개의 네 모칸이 나오는데 그것은 hys1, hys2, animal 각각 교차하여 차지하는 칸이다. 그런데 hys1은 첫째 형질의 모형에만 있으므로 R의 inverse 앞 뒤에 [ [ 1, 0], [0, 0] ]을 곱하여 준다. 구체적으로 다음과 같다.

같은 방식으로 RHS를 채우고, 나머지 개체들도 같은 방식으로 LHS와 RHS를 채운다.

# multiple trait animal model, unequal design matrices(different model) - Date : 2020-07-29

# Linear Models for the Prediction of Animal Breeding Values, 3rd Edtion

# Raphael Mrode

# Example 5.3 p80

위의 예를 가지고 R로 육종가를 구하고, 정확도를 구하고, 육종가를 분할하고, DYD를 계산한 것을 설명한다.

 

자료는 다음과 같다.

4 1 1 201 280
5 1 2 150 200
6 2 1 160 190
7 1 1 180 250
8 2 2 285 300

mt03_data.txt로 저장한다. 컬럼 설명은 위에서 하였다.

 

혈통은 다음과 같다.

1 0 0
2 0 0
3 0 0
4 1 2
5 3 2
6 1 5
7 3 4
8 1 7

mt3_pedi.txt로 저장한다.

 

R 프로그램은 다음과 같다.

# multiple trait animal model, unequal design matrices(different model) - Date : 2020-07-29

# Linear Models for the Prediction of Animal Breeding Values, 3rd Edtion

# Raphael Mrode

# Example 5.3 p80

# MASS packages
if (!("MASS" %in% installed.packages())) {
    install.packages("MASS", dependencies = TRUE)
}
library(MASS)

# functions

# find the position in mixed model equation(lhs and rhs)
pos_mme <- function(trts, lvls, vals) {
    pos = rep(0, length(vals))
    
    for (i in 1:length(vals)) {
        if (i == 1) {
            pos[i] = trts * (vals[i] - 1) + 1
        } else {
            pos[i] = trts * sum(lvls[1:i - 1]) + trts * (vals[i] - 1) + 1
        }
    }
    
    return(pos)
}

# 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(choose.dir())
setwd("D:/temp/06_multiple_traits_03_R")

# print working directory
getwd()

# no_of_trts
no_trts = 2

# list all possible combination of data
dtcombi0 = expand.grid(rep(list(0:1), no_trts))
dtcombi = dtcombi0[2:nrow(dtcombi0), ]
rownames(dtcombi) = NULL

# print
dtcombi

# prdigree and data

# read pedigree : animal sire dam
pedi = read.table("mt03_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("mt03_data.txt", header = FALSE, sep = "", col.names = c("animal", "hys1", "hys2", "fat1", "fat2"), 
    stringsAsFactors = FALSE)

# print
data

# number of data
no_data = nrow(data)

# print
no_data

# how many traits does animal have?
data2 = data.frame(data, dtsts = c(0))
data2$dtsts = ifelse(data2$fat1 != 0, data2$dtsts + 1, data2$dtsts)
data2$dtsts = ifelse(data2$fat2 != 0, data2$dtsts + 2, data2$dtsts)

# print
data2

# levels of animal, hys1, hys2
lvls_hys1 = max(data2$hys1)
lvls_hys2 = max(data2$hys2)
lvls_ani = max(data2$animal)

# print
lvls_hys1
lvls_hys2
lvls_ani

# effect data status. 1st effect(hys1) is for trait 1, 2nd effect(hys2) is for trait 2, 3rd effect is for trait 1 and 2.
effdtsts = array(rep(0, no_trts * no_trts * (2^no_trts - 1)), dim = c(no_trts, no_trts, (2^no_trts - 1)))
effdtsts[,,1] = diag(c(1,0))
effdtsts[,,2] = diag(c(0,1))
effdtsts[,,3] = diag(c(1,1))

# print
effdtsts

# variance component additive genetic
G = matrix(c(35, 28, 28, 30), 2, 2)
# residual
R = matrix(c(65, 27, 27, 70), 2, 2)

# print
G
R

# inverse of G
Gi = ginv(G)

# print
Gi

# inverse of R
Ri = array(rep(0, no_trts * no_trts * (2^no_trts - 1)), dim = c(no_trts, no_trts, (2^no_trts - 1)))

for (i in 1:(2^no_trts - 1)) {
    R0 = R
    R0[which(dtcombi[i, ] == 0), ] = 0
    R0[, which(dtcombi[i, ] == 0)] = 0
    Ri[, , i] = ginv(R0)
}

# print
Ri

# empty lhs
lhs = matrix(rep(0, (no_trts * (lvls_hys1 + lvls_hys2 + lvls_ani))^2), nrow = no_trts * (lvls_hys1 + lvls_hys2 + lvls_ani))

# print
dim(lhs)
lhs

# empty rhs
rhs = matrix(rep(0, (no_trts * (lvls_hys1 + lvls_hys2 + lvls_ani))), nrow = no_trts * (lvls_hys1 + lvls_hys2 + lvls_ani))

# print
dim(rhs)
rhs

# fill the MME
for (i in 1:no_data) {
#i = 1
    pos = pos_mme(no_trts, c(lvls_hys1, lvls_hys2, lvls_ani), c(data2$hys1[i], data2$hys2[i], data2$animal[i]))

    for (j in 1:length(pos)) {
        rfirst = pos[j]
        rlast = (pos[j] + no_trts - 1)

        for (k in 1:length(pos)) {
            cfirst = pos[k]
            clast = (pos[k] + no_trts - 1)
            
            lhs[rfirst : rlast, cfirst : clast] = lhs[rfirst : rlast, cfirst : clast] + effdtsts[,,j] %*% Ri[, , data2$dtsts[i]] %*% effdtsts[,,k]
        }
        rhs[rfirst : rlast] = rhs[rfirst : rlast] + effdtsts[,,j] %*% Ri[, , data2$dtsts[i]] %*% as.numeric(data2[i, 4:5])
    }
}

# print lhs and rhs
lhs
rhs

# inverse matrix of NRM
ai = ainv(pedi)

# print
ai

# add ai to lhs
afirst = no_trts * lvls_hys1 * lvls_hys2 + 1
alast = no_trts * (lvls_hys1 * lvls_hys2 + lvls_ani)
afirst
alast
lhs[afirst : alast, afirst : alast] = lhs[afirst : alast, afirst : alast] + ai %x% Gi


# print
#lhs[c(2,4,5,7),] = 0
#lhs[,c(2,4,5,7)] = 0
#rhs[c(2,4,5,7),] = 0

lhs
rhs

# generalised inverse of lhs
gi_lhs = ginv(lhs)

# print
gi_lhs

# solution
sol = gi_lhs %*% rhs

# print
sol

# levels of fixed effect 1 in traits 1
lvls_t1_f1 = rep(0, lvls_hys1)
for ( i in 1 : lvls_hys1){
    if (i == 1){
        lvls_t1_f1[i] = 1
    } else {
        lvls_t1_f1[i] = 1 + (i - 1) * no_trts
    }
}

# print
lvls_t1_f1

# levels of fixed effect 1 in traits 2
lvls_t2_f1 = lvls_t1_f1 + 1

# print
lvls_t2_f1

# levels of fixed effect 2 in traits 1
lvls_t1_f2 = rep(0, lvls_hys2)
for ( i in 1 : lvls_hys2){
    if (i == 1){
        lvls_t1_f2[i] = 1
    } else {
        lvls_t1_f2[i] = 1 + (i - 1) * no_trts
    }
}

# print
lvls_t1_f2 = lvls_t1_f2 + no_trts * lvls_hys1
lvls_t1_f2

# levels of fixed effect 2 in traits 2
lvls_t2_f2 = lvls_t1_f2 + 1

# print
lvls_t2_f2

# levels of animal effect in traits 1
lvls_t1_ani = rep(0, lvls_ani)
for ( i in 1 : lvls_ani){
    if (i == 1){
        lvls_t1_ani[i] = 1
    } else {
        lvls_t1_ani[i] = 1 + (i - 1) * no_trts
    }
}
lvls_t1_ani = lvls_t1_ani + no_trts * (lvls_hys1 + lvls_hys2)

# print
lvls_t1_ani

# levels of fixed effect in traits 1
lvls_t2_ani = lvls_t1_ani + 1

# print
lvls_t2_ani

# solutions for fixed effects 1
sol_t1_f1 = as.matrix(sol[lvls_t1_f1])
sol_t2_f1 = as.matrix(sol[lvls_t2_f1])
sol_f1 = cbind(sol_t1_f1, sol_t2_f1)

# print
sol_f1

# solutions for fixed effects 2
sol_t1_f2 = as.matrix(sol[lvls_t1_f2])
sol_t2_f2 = as.matrix(sol[lvls_t2_f2])
sol_f2 = cbind(sol_t1_f2, sol_t2_f2)

# print
sol_f2

# breedinv value
sol_t1_bv = sol[lvls_t1_ani]
sol_t2_bv = sol[lvls_t2_ani]
sol_bv = cbind(sol_t1_bv, sol_t2_bv)

# print
sol_bv

# reliability(r2), accuracy(r), standard error of prediction(SEP)

# diagonal elements of the generalized inverse of LHS for animal equation
d_ani_t1 = diag(gi_lhs[lvls_t1_ani, lvls_t1_ani])
d_ani_t2 = diag(gi_lhs[lvls_t2_ani, lvls_t2_ani])
d_ani = cbind(d_ani_t1, d_ani_t2)

# print
d_ani

# reliability
rel = matrix(rep(0, lvls_ani * no_trts), ncol = no_trts)

for (i in 1 : lvls_ani) {
    rel[i, ] = 1 - d_ani[i, ]/diag(G)
}

# print
rel

# accuracy
acc = sqrt(rel)

# print
acc

# standard error of prediction(SEP)
sep = sqrt(d_ani)

# 확인
sep

# partitioning of breeding values

# yield deviation
yd1 = matrix(rep(0, lvls_ani * no_trts), ncol = no_trts)

# numerator of n2
a2 = array(rep(0, no_trts * no_trts * lvls_ani), dim = c(no_trts, no_trts, lvls_ani))

# looping data
for (i in 1:no_data) {
    yd1[data[i, 1], ] = as.matrix(data2[i, c(4, 5)] - sol_f1[data2[i, 2], ] - sol_f2[data2[i, 3], ])
    
    a2[, , data[i, 1]] = Ri[, , data2$dtsts[i]]
}

# print
yd1
a2

# Parents average, progeny contribution

# parents avearge
pa1 = matrix(rep(0, lvls_ani * no_trts), ncol = no_trts)

# progeny contribution numerator
pc0 = matrix(rep(0, lvls_ani * no_trts), ncol = no_trts)

# numerator of n3, denominator of progeny contribution
a3 = array(rep(0, no_trts * no_trts * lvls_ani), dim = c(no_trts, no_trts, lvls_ani))

# numerator of n1
a1 = array(rep(0, no_trts * no_trts * lvls_ani), dim = c(no_trts, no_trts, lvls_ani))


# looping pedi
for (i in 1 : lvls_ani) {
    
    sire = pedi[i, 2]
    dam = pedi[i, 3]
    
    if (sire == 0 & dam == 0) {
        # both parents unknown PA
        a1[, , i] = 1 * Gi
        
    } else if (sire != 0 & dam == 0) {
        # 개체의 부만 알고 있을 경우
        
        # PA
        a1[, , i] = 4/3 * Gi
        pa1[i, ] = sol_bv[sire, ]/2
        
        # PC for sire
        a3[, , sire] = a3[, , sire] + 0.5 * Gi * (2/3)
        pc0[sire, ] = pc0[sire, ] + (0.5 * Gi * (2/3)) %*% (2 * sol_bv[i, ])
        
    } else if (sire == 0 & dam != 0) {
        # dam known
        
        # PA
        a1[, , i] = 4/3 * Gi
        pa1[i, ] = sol_bv[dam, ]/2
        
        # PC for dam
        a3[, , dam] = a3[, , dam] + 0.5 * Gi * (2/3)
        pc0[dam, ] = pc0[dam, ] + (0.5 * Gi * (2/3)) %*% (2 * sol_bv[i])
        
    } else {
        # both parents known
        
        # PA
        a1[, , i] = 2 * Gi
        pa1[i, ] = (sol_bv[sire, ] + sol_bv[dam, ])/2
        
        # PC for sire
        a3[, , sire] = a3[, , sire] + 0.5 * Gi
        pc0[sire, ] = pc0[sire, ] + (0.5 * Gi) %*% (2 * sol_bv[i, ] - sol_bv[dam, ])
        
        # PC for dam
        a3[, , dam] = a3[, , dam] + 0.5 * Gi
        pc0[dam, ] = pc0[dam, ] + (0.5 * Gi) %*% (2 * sol_bv[i, ] - sol_bv[sire, ])
        
    }
}

# print
a1
pa1
a3
pc0

# denominator of n1, n2, n3, diagonal of animals in LHS
n_de = a1 + a2 + a3

# print
n_de

# parents average fraction of breeding values
pa = matrix(rep(0, lvls_ani * no_trts), ncol = no_trts)
for (i in 1 : lvls_ani) {
    pa[i, ] = ginv(n_de[, , i]) %*% a1[, , i] %*% pa1[i, ]
}

# print
pa

# yield deviation fraction of breeding values
yd = matrix(rep(0, lvls_ani * no_trts), ncol = no_trts)
for (i in 1 : lvls_ani) {
    yd[i, ] = ginv(n_de[, , i]) %*% a2[, , i] %*% yd1[i, ]
}

# print
yd

# progeny contribution
pc1 = matrix(rep(0, lvls_ani * no_trts), ncol = no_trts)
for (i in 1 : lvls_ani) {
    pc1[i, ] = ginv(a3[, , i]) %*% pc0[i, ]
}
pc1[is.nan(pc1) == TRUE] = 0
pc1

# progeny contribution fraction of breeding values
pc = matrix(rep(0, lvls_ani * no_trts), ncol = no_trts)
for (i in 1 : lvls_ani) {
    pc[i, ] = ginv(n_de[, , i]) %*% a3[, , i] %*% pc1[i, ]
}

# print
pc

# Progeny(Daughter) Yield Deviation(PYD, DYD)

# n2 of progeny
n2prog = array(rep(0, no_trts * no_trts * lvls_ani), dim = c(no_trts, no_trts, lvls_ani))
for (i in 1 : lvls_ani) {
    n2prog[, , i] = ginv((a1 + a2)[, , i]) %*% a2[, , i]
}

# print
n2prog

# numerator of dyd : summation of u of progeny * n2 of progeny * (2 * YD - bv of mate)
dyd_n = matrix(rep(0, lvls_ani * no_trts), ncol = no_trts)
# denominator of dyd : summation of u of progeny * n2 of progeny
dyd_d = array(rep(0, no_trts * no_trts * lvls_ani), dim = c(no_trts, no_trts, lvls_ani))

# looping pedi
for (i in 1 : lvls_ani) {
    # i = 5
    sire = pedi[i, 2]
    dam = pedi[i, 3]
    
    if (sire == 0 & dam == 0) {
        # both parents unknown
        
    } else if (sire != 0 & dam == 0) {
        # 개체의 부만 알고 있을 경우
        
        # dyd_n
        dyd_n[sire, ] = dyd_n[sire, ] + (n2prog[, , i] * 2/3) %*% (2 * yd1[i, ])
        # dyd_d
        dyd_d[, , sire] = dyd_d[, , sire] + n2prog[, , i] * 2/3
        
    } else if (sire == 0 & dam != 0) {
        # dam known
        
        # dyd_n
        dyd_n[dam, ] = dyd_n[dam, ] + (n2prog[, , i] * 2/3) %*% (2 * yd1[i, ])
        # dyd_d
        dyd_d[, , dam] = dyd_d[, , dam] + n2prog[, , i] * 2/3
        
    } else {
        # both parents known
        
        # dyd_n
        dyd_n[sire, ] = dyd_n[sire, ] + n2prog[, , i] %*% (2 * yd1[i, ] - sol_bv[dam, ])
        dyd_n[dam, ] = dyd_n[dam, ] + n2prog[, , i] %*% (2 * yd1[i, ] - sol_bv[sire, ])
        
        # dyd_d
        dyd_d[, , sire] = dyd_d[, , sire] + n2prog[, , i]
        dyd_d[, , dam] = dyd_d[, , dam] + n2prog[, , i]
        
    }
}

# print
dyd_n
dyd_d

# dyd
dyd = matrix(rep(0, lvls_ani * no_trts), ncol = no_trts)

# looping pedi
for (i in 1 : lvls_ani) {
    dyd[i, ] = ginv(dyd_d[, , i]) %*% dyd_n[i, ]
}
dyd[is.nan(dyd) == TRUE] = 0

# print
dyd

# breeding values and fractions
mt3_result = data.frame(animal = pedi[, 1], animal_bv = sol_bv, rel = rel, acc = acc, sep = sep, pa = pa, 
    yd = yd, pc = pc, sum_of_fr = pa + yd + pc, dyd = dyd)

# print
mt3_result

# 파일로 출력, 분리자는 ',', 따옴표 없고
output_filename = gsub("[ -]", "", paste("mt3_result_", Sys.Date(), ".csv"))
write.table(mt3_result, output_filename, sep = ",", row.names = FALSE, quote = FALSE)

 

실행결과는 다음과 같다.

> # multiple trait animal model, unequal design matrices(different model) - Date : 2020-07-29
> 
> # Linear Models for the Prediction of Animal Breeding Values, 3rd Edtion
> 
> # Raphael Mrode
> 
> # Example 5.3 p80
> 
> # MASS packages
> if (!("MASS" %in% installed.packages())) {
+     install.packages("MASS", dependencies = TRUE)
+ }
> library(MASS)
> 
> # functions
> 
> # find the position in mixed model equation(lhs and rhs)
> pos_mme <- function(trts, lvls, vals) {
+     pos = rep(0, length(vals))
+     
+     for (i in 1:length(vals)) {
+         if (i == 1) {
+             pos[i] = trts * (vals[i] - 1) + 1
+         } else {
+             pos[i] = trts * sum(lvls[1:i - 1]) + trts * (vals[i] - 1) + 1
+         }
+     }
+     
+     return(pos)
+ }
> 
> # 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(choose.dir())
> setwd("D:/temp/06_multiple_traits_03_R")
> 
> # print working directory
> getwd()
[1] "D:/temp/06_multiple_traits_03_R"
> 
> # no_of_trts
> no_trts = 2
> 
> # list all possible combination of data
> dtcombi0 = expand.grid(rep(list(0:1), no_trts))
> dtcombi = dtcombi0[2:nrow(dtcombi0), ]
> rownames(dtcombi) = NULL
> 
> # print
> dtcombi
  Var1 Var2
1    1    0
2    0    1
3    1    1
> 
> # prdigree and data
> 
> # read pedigree : animal sire dam
> pedi = read.table("mt03_pedi.txt", header = FALSE, sep = "", col.names = c("animal", "sire", "dam"), 
+     stringsAsFactors = FALSE)
> pedi = pedi[order(pedi$animal), ]
> 
> # print
> pedi
  animal sire dam
1      1    0   0
2      2    0   0
3      3    0   0
4      4    1   2
5      5    3   2
6      6    1   5
7      7    3   4
8      8    1   7
> 
> # read data : animal, dam, sex, weaning_weight
> data = read.table("mt03_data.txt", header = FALSE, sep = "", col.names = c("animal", "hys1", "hys2", "fat1", "fat2"), 
+     stringsAsFactors = FALSE)
> 
> # print
> data
  animal hys1 hys2 fat1 fat2
1      4    1    1  201  280
2      5    1    2  150  200
3      6    2    1  160  190
4      7    1    1  180  250
5      8    2    2  285  300
> 
> # number of data
> no_data = nrow(data)
> 
> # print
> no_data
[1] 5
> 
> # how many traits does animal have?
> data2 = data.frame(data, dtsts = c(0))
> data2$dtsts = ifelse(data2$fat1 != 0, data2$dtsts + 1, data2$dtsts)
> data2$dtsts = ifelse(data2$fat2 != 0, data2$dtsts + 2, data2$dtsts)
> 
> # print
> data2
  animal hys1 hys2 fat1 fat2 dtsts
1      4    1    1  201  280     3
2      5    1    2  150  200     3
3      6    2    1  160  190     3
4      7    1    1  180  250     3
5      8    2    2  285  300     3
> 
> # levels of animal, hys1, hys2
> lvls_hys1 = max(data2$hys1)
> lvls_hys2 = max(data2$hys2)
> lvls_ani = max(data2$animal)
> 
> # print
> lvls_hys1
[1] 2
> lvls_hys2
[1] 2
> lvls_ani
[1] 8
> 
> # effect data status. 1st effect(hys1) is for trait 1, 2nd effect(hys2) is for trait 2, 3rd effect is for trait 1 and 2.
> effdtsts = array(rep(0, no_trts * no_trts * (2^no_trts - 1)), dim = c(no_trts, no_trts, (2^no_trts - 1)))
> effdtsts[,,1] = diag(c(1,0))
> effdtsts[,,2] = diag(c(0,1))
> effdtsts[,,3] = diag(c(1,1))
> 
> # print
> effdtsts
, , 1

     [,1] [,2]
[1,]    1    0
[2,]    0    0

, , 2

     [,1] [,2]
[1,]    0    0
[2,]    0    1

, , 3

     [,1] [,2]
[1,]    1    0
[2,]    0    1

> 
> # variance component additive genetic
> G = matrix(c(35, 28, 28, 30), 2, 2)
> # residual
> R = matrix(c(65, 27, 27, 70), 2, 2)
> 
> # print
> G
     [,1] [,2]
[1,]   35   28
[2,]   28   30
> R
     [,1] [,2]
[1,]   65   27
[2,]   27   70
> 
> # inverse of G
> Gi = ginv(G)
> 
> # print
> Gi
           [,1]       [,2]
[1,]  0.1127820 -0.1052632
[2,] -0.1052632  0.1315789
> 
> # inverse of R
> Ri = array(rep(0, no_trts * no_trts * (2^no_trts - 1)), dim = c(no_trts, no_trts, (2^no_trts - 1)))
> 
> for (i in 1:(2^no_trts - 1)) {
+     R0 = R
+     R0[which(dtcombi[i, ] == 0), ] = 0
+     R0[, which(dtcombi[i, ] == 0)] = 0
+     Ri[, , i] = ginv(R0)
+ }
> 
> # print
> Ri
, , 1

           [,1] [,2]
[1,] 0.01538462    0
[2,] 0.00000000    0

, , 2

     [,1]       [,2]
[1,]    0 0.00000000
[2,]    0 0.01428571

, , 3

             [,1]         [,2]
[1,]  0.018319812 -0.007066213
[2,] -0.007066213  0.017011254

> 
> # empty lhs
> lhs = matrix(rep(0, (no_trts * (lvls_hys1 + lvls_hys2 + lvls_ani))^2), nrow = no_trts * (lvls_hys1 + lvls_hys2 + lvls_ani))
> 
> # print
> dim(lhs)
[1] 24 24
> lhs
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21]
 [1,]    0    0    0    0    0    0    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     0     0     0     0     0     0
 [3,]    0    0    0    0    0    0    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     0     0     0     0     0     0
 [5,]    0    0    0    0    0    0    0    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     0     0     0     0     0     0     0
 [7,]    0    0    0    0    0    0    0    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     0     0     0     0     0     0     0
 [9,]    0    0    0    0    0    0    0    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     0     0     0     0     0     0     0
[11,]    0    0    0    0    0    0    0    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     0     0     0     0     0     0     0
[13,]    0    0    0    0    0    0    0    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     0     0     0     0     0     0     0
[15,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0     0
[16,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0     0
[17,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0     0
[18,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0     0
[19,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0     0
[20,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0     0
[21,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0     0
[22,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0     0
[23,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0     0
[24,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0     0
      [,22] [,23] [,24]
 [1,]     0     0     0
 [2,]     0     0     0
 [3,]     0     0     0
 [4,]     0     0     0
 [5,]     0     0     0
 [6,]     0     0     0
 [7,]     0     0     0
 [8,]     0     0     0
 [9,]     0     0     0
[10,]     0     0     0
[11,]     0     0     0
[12,]     0     0     0
[13,]     0     0     0
[14,]     0     0     0
[15,]     0     0     0
[16,]     0     0     0
[17,]     0     0     0
[18,]     0     0     0
[19,]     0     0     0
[20,]     0     0     0
[21,]     0     0     0
[22,]     0     0     0
[23,]     0     0     0
[24,]     0     0     0
> 
> # empty rhs
> rhs = matrix(rep(0, (no_trts * (lvls_hys1 + lvls_hys2 + lvls_ani))), nrow = no_trts * (lvls_hys1 + lvls_hys2 + lvls_ani))
> 
> # print
> dim(rhs)
[1] 24  1
> 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
[15,]    0
[16,]    0
[17,]    0
[18,]    0
[19,]    0
[20,]    0
[21,]    0
[22,]    0
[23,]    0
[24,]    0
> 
> # fill the MME
> for (i in 1:no_data) {
+ #i = 1
+     pos = pos_mme(no_trts, c(lvls_hys1, lvls_hys2, lvls_ani), c(data2$hys1[i], data2$hys2[i], data2$animal[i]))
+ 
+     for (j in 1:length(pos)) {
+         rfirst = pos[j]
+         rlast = (pos[j] + no_trts - 1)
+ 
+         for (k in 1:length(pos)) {
+             cfirst = pos[k]
+             clast = (pos[k] + no_trts - 1)
+             
+             lhs[rfirst : rlast, cfirst : clast] = lhs[rfirst : rlast, cfirst : clast] + effdtsts[,,j] %*% Ri[, , data2$dtsts[i]] %*% effdtsts[,,k]
+         }
+         rhs[rfirst : rlast] = rhs[rfirst : rlast] + effdtsts[,,j] %*% Ri[, , data2$dtsts[i]] %*% as.numeric(data2[i, 4:5])
+     }
+ }
> 
> # print lhs and rhs
> lhs
              [,1] [,2]         [,3] [,4] [,5]         [,6] [,7]         [,8] [,9] [,10] [,11] [,12] [,13] [,14]
 [1,]  0.054959435    0  0.000000000    0    0 -0.014132426    0 -0.007066213    0     0     0     0     0     0
 [2,]  0.000000000    0  0.000000000    0    0  0.000000000    0  0.000000000    0     0     0     0     0     0
 [3,]  0.000000000    0  0.036639623    0    0 -0.007066213    0 -0.007066213    0     0     0     0     0     0
 [4,]  0.000000000    0  0.000000000    0    0  0.000000000    0  0.000000000    0     0     0     0     0     0
 [5,]  0.000000000    0  0.000000000    0    0  0.000000000    0  0.000000000    0     0     0     0     0     0
 [6,] -0.014132426    0 -0.007066213    0    0  0.051033761    0  0.000000000    0     0     0     0     0     0
 [7,]  0.000000000    0  0.000000000    0    0  0.000000000    0  0.000000000    0     0     0     0     0     0
 [8,] -0.007066213    0 -0.007066213    0    0  0.000000000    0  0.034022507    0     0     0     0     0     0
 [9,]  0.000000000    0  0.000000000    0    0  0.000000000    0  0.000000000    0     0     0     0     0     0
[10,]  0.000000000    0  0.000000000    0    0  0.000000000    0  0.000000000    0     0     0     0     0     0
[11,]  0.000000000    0  0.000000000    0    0  0.000000000    0  0.000000000    0     0     0     0     0     0
[12,]  0.000000000    0  0.000000000    0    0  0.000000000    0  0.000000000    0     0     0     0     0     0
[13,]  0.000000000    0  0.000000000    0    0  0.000000000    0  0.000000000    0     0     0     0     0     0
[14,]  0.000000000    0  0.000000000    0    0  0.000000000    0  0.000000000    0     0     0     0     0     0
[15,]  0.018319812    0  0.000000000    0    0 -0.007066213    0  0.000000000    0     0     0     0     0     0
[16,] -0.007066213    0  0.000000000    0    0  0.017011254    0  0.000000000    0     0     0     0     0     0
[17,]  0.018319812    0  0.000000000    0    0  0.000000000    0 -0.007066213    0     0     0     0     0     0
[18,] -0.007066213    0  0.000000000    0    0  0.000000000    0  0.017011254    0     0     0     0     0     0
[19,]  0.000000000    0  0.018319812    0    0 -0.007066213    0  0.000000000    0     0     0     0     0     0
[20,]  0.000000000    0 -0.007066213    0    0  0.017011254    0  0.000000000    0     0     0     0     0     0
[21,]  0.018319812    0  0.000000000    0    0 -0.007066213    0  0.000000000    0     0     0     0     0     0
[22,] -0.007066213    0  0.000000000    0    0  0.017011254    0  0.000000000    0     0     0     0     0     0
[23,]  0.000000000    0  0.018319812    0    0  0.000000000    0 -0.007066213    0     0     0     0     0     0
[24,]  0.000000000    0 -0.007066213    0    0  0.000000000    0  0.017011254    0     0     0     0     0     0
             [,15]        [,16]        [,17]        [,18]        [,19]        [,20]        [,21]        [,22]        [,23]
 [1,]  0.018319812 -0.007066213  0.018319812 -0.007066213  0.000000000  0.000000000  0.018319812 -0.007066213  0.000000000
 [2,]  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000
 [3,]  0.000000000  0.000000000  0.000000000  0.000000000  0.018319812 -0.007066213  0.000000000  0.000000000  0.018319812
 [4,]  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000
 [5,]  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000
 [6,] -0.007066213  0.017011254  0.000000000  0.000000000 -0.007066213  0.017011254 -0.007066213  0.017011254  0.000000000
 [7,]  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000
 [8,]  0.000000000  0.000000000 -0.007066213  0.017011254  0.000000000  0.000000000  0.000000000  0.000000000 -0.007066213
 [9,]  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000
[10,]  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000
[11,]  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000
[12,]  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000
[13,]  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000
[14,]  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000
[15,]  0.018319812 -0.007066213  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000
[16,] -0.007066213  0.017011254  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000
[17,]  0.000000000  0.000000000  0.018319812 -0.007066213  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000
[18,]  0.000000000  0.000000000 -0.007066213  0.017011254  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000
[19,]  0.000000000  0.000000000  0.000000000  0.000000000  0.018319812 -0.007066213  0.000000000  0.000000000  0.000000000
[20,]  0.000000000  0.000000000  0.000000000  0.000000000 -0.007066213  0.017011254  0.000000000  0.000000000  0.000000000
[21,]  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.018319812 -0.007066213  0.000000000
[22,]  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000 -0.007066213  0.017011254  0.000000000
[23,]  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.018319812
[24,]  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000 -0.007066213
             [,24]
 [1,]  0.000000000
 [2,]  0.000000000
 [3,] -0.007066213
 [4,]  0.000000000
 [5,]  0.000000000
 [6,]  0.000000000
 [7,]  0.000000000
 [8,]  0.017011254
 [9,]  0.000000000
[10,]  0.000000000
[11,]  0.000000000
[12,]  0.000000000
[13,]  0.000000000
[14,]  0.000000000
[15,]  0.000000000
[16,]  0.000000000
[17,]  0.000000000
[18,]  0.000000000
[19,]  0.000000000
[20,]  0.000000000
[21,]  0.000000000
[22,]  0.000000000
[23,] -0.007066213
[24,]  0.017011254
> rhs
          [,1]
 [1,] 4.569484
 [2,] 0.000000
 [3,] 4.689872
 [4,] 0.000000
 [5,] 0.000000
 [6,] 8.425281
 [7,] 0.000000
 [8,] 5.431824
 [9,] 0.000000
[10,] 0.000000
[11,] 0.000000
[12,] 0.000000
[13,] 0.000000
[14,] 0.000000
[15,] 1.703742
[16,] 3.342842
[17,] 1.334729
[18,] 2.342319
[19,] 1.588589
[20,] 2.101544
[21,] 1.531013
[22,] 2.980895
[23,] 3.101282
[24,] 3.089505
> 
> # inverse matrix of NRM
> ai = ainv(pedi)
> 
> # print
> ai
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]  2.5  0.5  0.0 -1.0  0.5   -1  0.5   -1
[2,]  0.5  2.0  0.5 -1.0 -1.0    0  0.0    0
[3,]  0.0  0.5  2.0  0.5 -1.0    0 -1.0    0
[4,] -1.0 -1.0  0.5  2.5  0.0    0 -1.0    0
[5,]  0.5 -1.0 -1.0  0.0  2.5   -1  0.0    0
[6,] -1.0  0.0  0.0  0.0 -1.0    2  0.0    0
[7,]  0.5  0.0 -1.0 -1.0  0.0    0  2.5   -1
[8,] -1.0  0.0  0.0  0.0  0.0    0 -1.0    2
> 
> # add ai to lhs
> afirst = no_trts * lvls_hys1 * lvls_hys2 + 1
> alast = no_trts * (lvls_hys1 * lvls_hys2 + lvls_ani)
> afirst
[1] 9
> alast
[1] 24
> lhs[afirst : alast, afirst : alast] = lhs[afirst : alast, afirst : alast] + ai %x% Gi
> 
> 
> # print
> #lhs[c(2,4,5,7),] = 0
> #lhs[,c(2,4,5,7)] = 0
> #rhs[c(2,4,5,7),] = 0
> 
> lhs
              [,1] [,2]         [,3] [,4] [,5]         [,6] [,7]         [,8]        [,9]       [,10]       [,11]
 [1,]  0.054959435    0  0.000000000    0    0 -0.014132426    0 -0.007066213  0.00000000  0.00000000  0.00000000
 [2,]  0.000000000    0  0.000000000    0    0  0.000000000    0  0.000000000  0.00000000  0.00000000  0.00000000
 [3,]  0.000000000    0  0.036639623    0    0 -0.007066213    0 -0.007066213  0.00000000  0.00000000  0.00000000
 [4,]  0.000000000    0  0.000000000    0    0  0.000000000    0  0.000000000  0.00000000  0.00000000  0.00000000
 [5,]  0.000000000    0  0.000000000    0    0  0.000000000    0  0.000000000  0.00000000  0.00000000  0.00000000
 [6,] -0.014132426    0 -0.007066213    0    0  0.051033761    0  0.000000000  0.00000000  0.00000000  0.00000000
 [7,]  0.000000000    0  0.000000000    0    0  0.000000000    0  0.000000000  0.00000000  0.00000000  0.00000000
 [8,] -0.007066213    0 -0.007066213    0    0  0.000000000    0  0.034022507  0.00000000  0.00000000  0.00000000
 [9,]  0.000000000    0  0.000000000    0    0  0.000000000    0  0.000000000  0.28195489 -0.26315789  0.05639098
[10,]  0.000000000    0  0.000000000    0    0  0.000000000    0  0.000000000 -0.26315789  0.32894737 -0.05263158
[11,]  0.000000000    0  0.000000000    0    0  0.000000000    0  0.000000000  0.05639098 -0.05263158  0.22556391
[12,]  0.000000000    0  0.000000000    0    0  0.000000000    0  0.000000000 -0.05263158  0.06578947 -0.21052632
[13,]  0.000000000    0  0.000000000    0    0  0.000000000    0  0.000000000  0.00000000  0.00000000  0.05639098
[14,]  0.000000000    0  0.000000000    0    0  0.000000000    0  0.000000000  0.00000000  0.00000000 -0.05263158
[15,]  0.018319812    0  0.000000000    0    0 -0.007066213    0  0.000000000 -0.11278195  0.10526316 -0.11278195
[16,] -0.007066213    0  0.000000000    0    0  0.017011254    0  0.000000000  0.10526316 -0.13157895  0.10526316
[17,]  0.018319812    0  0.000000000    0    0  0.000000000    0 -0.007066213  0.05639098 -0.05263158 -0.11278195
[18,] -0.007066213    0  0.000000000    0    0  0.000000000    0  0.017011254 -0.05263158  0.06578947  0.10526316
[19,]  0.000000000    0  0.018319812    0    0 -0.007066213    0  0.000000000 -0.11278195  0.10526316  0.00000000
[20,]  0.000000000    0 -0.007066213    0    0  0.017011254    0  0.000000000  0.10526316 -0.13157895  0.00000000
[21,]  0.018319812    0  0.000000000    0    0 -0.007066213    0  0.000000000  0.05639098 -0.05263158  0.00000000
[22,] -0.007066213    0  0.000000000    0    0  0.017011254    0  0.000000000 -0.05263158  0.06578947  0.00000000
[23,]  0.000000000    0  0.018319812    0    0  0.000000000    0 -0.007066213 -0.11278195  0.10526316  0.00000000
[24,]  0.000000000    0 -0.007066213    0    0  0.000000000    0  0.017011254  0.10526316 -0.13157895  0.00000000
            [,12]       [,13]       [,14]        [,15]        [,16]        [,17]        [,18]        [,19]        [,20]
 [1,]  0.00000000  0.00000000  0.00000000  0.018319812 -0.007066213  0.018319812 -0.007066213  0.000000000  0.000000000
 [2,]  0.00000000  0.00000000  0.00000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000
 [3,]  0.00000000  0.00000000  0.00000000  0.000000000  0.000000000  0.000000000  0.000000000  0.018319812 -0.007066213
 [4,]  0.00000000  0.00000000  0.00000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000
 [5,]  0.00000000  0.00000000  0.00000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000
 [6,]  0.00000000  0.00000000  0.00000000 -0.007066213  0.017011254  0.000000000  0.000000000 -0.007066213  0.017011254
 [7,]  0.00000000  0.00000000  0.00000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000
 [8,]  0.00000000  0.00000000  0.00000000  0.000000000  0.000000000 -0.007066213  0.017011254  0.000000000  0.000000000
 [9,] -0.05263158  0.00000000  0.00000000 -0.112781955  0.105263158  0.056390977 -0.052631579 -0.112781955  0.105263158
[10,]  0.06578947  0.00000000  0.00000000  0.105263158 -0.131578947 -0.052631579  0.065789474  0.105263158 -0.131578947
[11,] -0.21052632  0.05639098 -0.05263158 -0.112781955  0.105263158 -0.112781955  0.105263158  0.000000000  0.000000000
[12,]  0.26315789 -0.05263158  0.06578947  0.105263158 -0.131578947  0.105263158 -0.131578947  0.000000000  0.000000000
[13,] -0.05263158  0.22556391 -0.21052632  0.056390977 -0.052631579 -0.112781955  0.105263158  0.000000000  0.000000000
[14,]  0.06578947 -0.21052632  0.26315789 -0.052631579  0.065789474  0.105263158 -0.131578947  0.000000000  0.000000000
[15,]  0.10526316  0.05639098 -0.05263158  0.300274699 -0.270224108  0.000000000  0.000000000  0.000000000  0.000000000
[16,] -0.13157895 -0.05263158  0.06578947 -0.270224108  0.345958622  0.000000000  0.000000000  0.000000000  0.000000000
[17,]  0.10526316 -0.11278195  0.10526316  0.000000000  0.000000000  0.300274699 -0.270224108 -0.112781955  0.105263158
[18,] -0.13157895  0.10526316 -0.13157895  0.000000000  0.000000000 -0.270224108  0.345958622  0.105263158 -0.131578947
[19,]  0.00000000  0.00000000  0.00000000  0.000000000  0.000000000 -0.112781955  0.105263158  0.243883721 -0.217592529
[20,]  0.00000000  0.00000000  0.00000000  0.000000000  0.000000000  0.105263158 -0.131578947 -0.217592529  0.280169148
[21,]  0.00000000 -0.11278195  0.10526316 -0.112781955  0.105263158  0.000000000  0.000000000  0.000000000  0.000000000
[22,]  0.00000000  0.10526316 -0.13157895  0.105263158 -0.131578947  0.000000000  0.000000000  0.000000000  0.000000000
[23,]  0.00000000  0.00000000  0.00000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000
[24,]  0.00000000  0.00000000  0.00000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000
             [,21]        [,22]        [,23]        [,24]
 [1,]  0.018319812 -0.007066213  0.000000000  0.000000000
 [2,]  0.000000000  0.000000000  0.000000000  0.000000000
 [3,]  0.000000000  0.000000000  0.018319812 -0.007066213
 [4,]  0.000000000  0.000000000  0.000000000  0.000000000
 [5,]  0.000000000  0.000000000  0.000000000  0.000000000
 [6,] -0.007066213  0.017011254  0.000000000  0.000000000
 [7,]  0.000000000  0.000000000  0.000000000  0.000000000
 [8,]  0.000000000  0.000000000 -0.007066213  0.017011254
 [9,]  0.056390977 -0.052631579 -0.112781955  0.105263158
[10,] -0.052631579  0.065789474  0.105263158 -0.131578947
[11,]  0.000000000  0.000000000  0.000000000  0.000000000
[12,]  0.000000000  0.000000000  0.000000000  0.000000000
[13,] -0.112781955  0.105263158  0.000000000  0.000000000
[14,]  0.105263158 -0.131578947  0.000000000  0.000000000
[15,] -0.112781955  0.105263158  0.000000000  0.000000000
[16,]  0.105263158 -0.131578947  0.000000000  0.000000000
[17,]  0.000000000  0.000000000  0.000000000  0.000000000
[18,]  0.000000000  0.000000000  0.000000000  0.000000000
[19,]  0.000000000  0.000000000  0.000000000  0.000000000
[20,]  0.000000000  0.000000000  0.000000000  0.000000000
[21,]  0.300274699 -0.270224108 -0.112781955  0.105263158
[22,] -0.270224108  0.345958622  0.105263158 -0.131578947
[23,] -0.112781955  0.105263158  0.243883721 -0.217592529
[24,]  0.105263158 -0.131578947 -0.217592529  0.280169148
> rhs
          [,1]
 [1,] 4.569484
 [2,] 0.000000
 [3,] 4.689872
 [4,] 0.000000
 [5,] 0.000000
 [6,] 8.425281
 [7,] 0.000000
 [8,] 5.431824
 [9,] 0.000000
[10,] 0.000000
[11,] 0.000000
[12,] 0.000000
[13,] 0.000000
[14,] 0.000000
[15,] 1.703742
[16,] 3.342842
[17,] 1.334729
[18,] 2.342319
[19,] 1.588589
[20,] 2.101544
[21,] 1.531013
[22,] 2.980895
[23,] 3.101282
[24,] 3.089505
> 
> # generalised inverse of lhs
> gi_lhs = ginv(lhs)
> 
> # print
> gi_lhs
               [,1]          [,2]          [,3]          [,4]          [,5]          [,6]          [,7]          [,8]
 [1,]  3.951667e+01 -1.239267e-12  1.860604e+01 -6.509892e-13  3.799885e-14  2.072396e+01  2.796365e-13  1.881434e+01
 [2,]  3.479718e-13 -2.059779e-26  4.537402e-13 -1.100182e-26  4.135091e-28  3.395211e-13  4.885042e-27  3.596606e-13
 [3,]  1.860604e+01 -1.245370e-12  5.131448e+01 -6.699507e-13  3.816029e-14  1.929738e+01  3.070411e-13  2.174017e+01
 [4,]  5.008545e-13 -2.928304e-26  5.333186e-13 -1.534105e-26  9.697558e-28  4.650294e-13  6.814292e-27  4.751315e-13
 [5,] -4.191738e-14  4.322369e-28  4.630835e-14  2.205518e-28  1.635113e-29 -2.098439e-14 -2.103332e-29  7.471397e-15
 [6,]  2.072396e+01 -1.143524e-12  1.929738e+01 -6.076725e-13  3.531775e-14  3.850489e+01  2.621040e-13  1.690015e+01
 [7,]  1.221833e-13 -5.346853e-27  9.330723e-14 -2.841939e-27  1.145413e-28  9.919806e-14  1.213193e-27  9.155333e-14
 [8,]  1.881434e+01 -1.174494e-12  2.174017e+01 -6.268728e-13  4.638977e-14  1.690015e+01  2.754894e-13  4.857570e+01
 [9,] -1.021391e+01  7.449606e-13 -1.667356e+01  4.070415e-13 -2.923215e-15 -1.037145e+01 -1.909223e-13 -9.939212e+00
[10,] -8.632111e+00  6.800524e-13 -1.255250e+01  3.699306e-13 -3.392157e-15 -1.143780e+01 -1.734011e-13 -1.024113e+01
[11,] -1.359657e+01  7.337009e-13 -8.734513e+00  3.996451e-13 -1.174280e-15 -9.572888e+00 -1.597105e-13 -9.030811e+00
[12,] -1.056770e+01  6.674130e-13 -7.534315e+00  3.637596e-13 -2.417536e-15 -1.020152e+01 -1.471988e-13 -9.694291e+00
[13,] -1.118952e+01  7.488013e-13 -9.591926e+00  3.768026e-13 -8.005536e-14 -8.055662e+00 -1.674551e-13 -9.029977e+00
[14,] -8.800189e+00  6.834317e-13 -7.913183e+00  3.431611e-13 -7.547304e-14 -8.360686e+00 -1.524225e-13 -1.006458e+01
[15,] -1.955071e+01  1.061556e-12 -1.609236e+01  5.613977e-13 -4.101809e-15 -1.560657e+01 -2.418968e-13 -1.292920e+01
[16,] -1.548013e+01  9.673848e-13 -1.311173e+01  5.112720e-13 -5.505309e-15 -1.718204e+01 -2.212832e-13 -1.324989e+01
[17,] -1.834414e+01  1.026361e-12 -1.450941e+01  5.449074e-13 -5.382512e-14 -1.275277e+01 -2.311228e-13 -1.461702e+01
[18,] -1.437142e+01  9.366364e-13 -1.218974e+01  4.965954e-13 -5.183407e-14 -1.312024e+01 -2.119496e-13 -1.629155e+01
[19,] -1.577718e+01  1.046802e-12 -2.228063e+01  5.707973e-13 -3.561339e-14 -1.503955e+01 -2.553750e-13 -1.419753e+01
[20,] -1.311613e+01  9.546899e-13 -1.708580e+01  5.190311e-13 -3.423878e-14 -1.637705e+01 -2.322789e-13 -1.473904e+01
[21,] -2.060855e+01  1.174903e-12 -1.708788e+01  5.898121e-13 -4.790546e-14 -1.594829e+01 -2.629515e-13 -1.442294e+01
[22,] -1.625287e+01  1.071616e-12 -1.395965e+01  5.369640e-13 -4.713575e-14 -1.729292e+01 -2.398455e-13 -1.530969e+01
[23,] -1.648150e+01  1.117432e-12 -2.347681e+01  5.676807e-13 -2.817260e-14 -1.441947e+01 -2.663465e-13 -1.675668e+01
[24,] -1.358002e+01  1.021215e-12 -1.802558e+01  5.161925e-13 -2.867988e-14 -1.534274e+01 -2.421319e-13 -1.826167e+01
               [,9]         [,10]         [,11]         [,12]         [,13]         [,14]         [,15]         [,16]
 [1,] -1.021391e+01 -8.632111e+00 -1.359657e+01 -1.056770e+01 -1.118952e+01 -8.800189e+00 -1.955071e+01 -1.548013e+01
 [2,] -2.934061e-13 -2.609547e-13 -2.427449e-13 -2.251661e-13 -9.650313e-14 -8.767265e-14 -3.475581e-13 -3.168947e-13
 [3,] -1.667356e+01 -1.255250e+01 -8.734513e+00 -7.534315e+00 -9.591926e+00 -7.913183e+00 -1.609236e+01 -1.311173e+01
 [4,] -3.433378e-13 -3.114666e-13 -2.871248e-13 -2.598961e-13 -2.856351e-13 -2.586731e-13 -4.366922e-13 -3.962122e-13
 [5,] -3.467913e-14 -2.779537e-14  5.074145e-14  4.160815e-14 -6.537701e-16 -4.172073e-15  3.000366e-14  2.502117e-14
 [6,] -1.037145e+01 -1.143780e+01 -9.572888e+00 -1.020152e+01 -8.055662e+00 -8.360686e+00 -1.560657e+01 -1.718204e+01
 [7,] -3.750542e-14 -3.166236e-14 -8.643645e-14 -7.331050e-14 -2.616409e-14 -1.882147e-14 -9.517557e-14 -8.096461e-14
 [8,] -9.939212e+00 -1.024113e+01 -9.030811e+00 -9.694291e+00 -9.029977e+00 -1.006458e+01 -1.292920e+01 -1.324989e+01
 [9,]  3.209312e+01  2.526634e+01  7.754462e-01  8.251871e-01  2.131438e+00  1.908472e+00  1.518990e+01  1.209336e+01
[10,]  2.526634e+01  2.731163e+01  8.346382e-01  8.935145e-01  1.899021e+00  1.794857e+00  1.212683e+01  1.335349e+01
[11,]  7.754462e-01  8.346382e-01  3.368775e+01  2.674286e+01  5.368049e-01  4.224971e-01  1.680003e+01  1.333016e+01
[12,]  8.251871e-01  8.935145e-01  2.674286e+01  2.874419e+01  4.319481e-01  3.622908e-01  1.331245e+01  1.432298e+01
[13,]  2.131438e+00  1.899021e+00  5.368049e-01  4.319481e-01  3.233176e+01  2.566903e+01  3.010069e+00  2.576479e+00
[14,]  1.908472e+00  1.794857e+00  4.224971e-01  3.622908e-01  2.566903e+01  2.784285e+01  2.560717e+00  2.323531e+00
[15,]  1.518990e+01  1.212683e+01  1.680003e+01  1.331245e+01  3.010069e+00  2.560717e+00  3.044248e+01  2.405106e+01
[16,]  1.209336e+01  1.335349e+01  1.333016e+01  1.432298e+01  2.576479e+00  2.323531e+00  2.405106e+01  2.638380e+01
[17,]  3.473266e+00  3.125123e+00  1.623159e+01  1.280185e+01  1.529514e+01  1.207303e+01  1.225756e+01  9.944174e+00
[18,]  3.144419e+00  2.986782e+00  1.278414e+01  1.379331e+01  1.207144e+01  1.321991e+01  9.917610e+00  1.010066e+01
[19,]  1.727846e+01  1.360962e+01  8.769870e+00  7.171554e+00  8.951673e+00  7.218822e+00  1.450576e+01  1.182501e+01
[20,]  1.360858e+01  1.451010e+01  7.178216e+00  7.776840e+00  7.213206e+00  7.713058e+00  1.184820e+01  1.263167e+01
[21,]  8.772284e+00  7.153655e+00  1.008590e+01  8.089705e+00  1.614181e+01  1.275664e+01  1.738384e+01  1.389940e+01
[22,]  7.131799e+00  7.726434e+00  8.097368e+00  8.464855e+00  1.277083e+01  1.380871e+01  1.388558e+01  1.489979e+01
[23,]  1.927493e+01  1.498610e+01  6.371341e+00  5.396174e+00  9.353731e+00  7.617724e+00  1.624722e+01  1.294420e+01
[24,]  1.501462e+01  1.621884e+01  5.390406e+00  5.635149e+00  7.594978e+00  8.146006e+00  1.297565e+01  1.402997e+01
              [,17]         [,18]         [,19]         [,20]         [,21]         [,22]         [,23]         [,24]
 [1,] -1.834414e+01 -1.437142e+01 -1.577718e+01 -1.311613e+01 -2.060855e+01 -1.625287e+01 -1.648150e+01 -1.358002e+01
 [2,] -2.558406e-13 -2.351251e-13 -3.425910e-13 -3.061498e-13 -2.813036e-13 -2.556802e-13 -3.543026e-13 -3.162754e-13
 [3,] -1.450941e+01 -1.218974e+01 -2.228063e+01 -1.708580e+01 -1.708788e+01 -1.395965e+01 -2.347681e+01 -1.802558e+01
 [4,] -3.945893e-13 -3.573265e-13 -4.372485e-13 -3.954858e-13 -4.739696e-13 -4.298512e-13 -4.908237e-13 -4.446901e-13
 [5,]  2.548488e-14  1.788768e-14 -1.943013e-14 -1.513501e-14  1.781725e-14  1.244479e-14 -2.986322e-14 -2.546035e-14
 [6,] -1.275277e+01 -1.312024e+01 -1.503955e+01 -1.637705e+01 -1.594829e+01 -1.729292e+01 -1.441947e+01 -1.534274e+01
 [7,] -8.435560e-14 -6.964001e-14 -7.292402e-14 -6.163944e-14 -8.394071e-14 -6.995659e-14 -6.886013e-14 -5.828654e-14
 [8,] -1.461702e+01 -1.629155e+01 -1.419753e+01 -1.473904e+01 -1.442294e+01 -1.530969e+01 -1.675668e+01 -1.826167e+01
 [9,]  3.473266e+00  3.144419e+00  1.727846e+01  1.360858e+01  8.772284e+00  7.131799e+00  1.927493e+01  1.501462e+01
[10,]  3.125123e+00  2.986782e+00  1.360962e+01  1.451010e+01  7.153655e+00  7.726434e+00  1.498610e+01  1.621884e+01
[11,]  1.623159e+01  1.278414e+01  8.769870e+00  7.178216e+00  1.008590e+01  8.097368e+00  6.371341e+00  5.390406e+00
[12,]  1.280185e+01  1.379331e+01  7.171554e+00  7.776840e+00  8.089705e+00  8.464855e+00  5.396174e+00  5.635149e+00
[13,]  1.529514e+01  1.207144e+01  8.951673e+00  7.213206e+00  1.614181e+01  1.277083e+01  9.353731e+00  7.594978e+00
[14,]  1.207303e+01  1.321991e+01  7.218822e+00  7.713058e+00  1.275664e+01  1.380871e+01  7.617724e+00  8.146006e+00
[15,]  1.225756e+01  9.917610e+00  1.450576e+01  1.184820e+01  1.738384e+01  1.388558e+01  1.624722e+01  1.297565e+01
[16,]  9.944174e+00  1.010066e+01  1.182501e+01  1.263167e+01  1.389940e+01  1.489979e+01  1.294420e+01  1.402997e+01
[17,]  2.958983e+01  2.325860e+01  1.614904e+01  1.291913e+01  1.524502e+01  1.226047e+01  1.080979e+01  9.109963e+00
[18,]  2.325860e+01  2.558931e+01  1.293232e+01  1.403359e+01  1.223516e+01  1.279749e+01  9.150064e+00  9.422754e+00
[19,]  1.614904e+01  1.293232e+01  3.116542e+01  2.446180e+01  1.339212e+01  1.100360e+01  1.668047e+01  1.329099e+01
[20,]  1.291913e+01  1.403359e+01  2.446180e+01  2.656940e+01  1.102049e+01  1.159678e+01  1.327036e+01  1.377779e+01
[21,]  1.524502e+01  1.223516e+01  1.339212e+01  1.102049e+01  3.077347e+01  2.427267e+01  1.920696e+01  1.526302e+01
[22,]  1.226047e+01  1.279749e+01  1.100360e+01  1.159678e+01  2.427267e+01  2.650225e+01  1.525557e+01  1.670182e+01
[23,]  1.080979e+01  9.150064e+00  1.668047e+01  1.327036e+01  1.920696e+01  1.525557e+01  3.345367e+01  2.615160e+01
[24,]  9.109963e+00  9.422754e+00  1.329099e+01  1.377779e+01  1.526302e+01  1.670182e+01  2.615160e+01  2.861921e+01
> 
> # solution
> sol = gi_lhs %*% rhs
> 
> # print
> sol
               [,1]
 [1,]  1.757313e+02
 [2,]  1.532129e-12
 [3,]  2.196133e+02
 [4,]  1.427759e-12
 [5,] -6.947473e-14
 [6,]  2.432391e+02
 [7,]  6.444061e-13
 [8,]  2.405497e+02
 [9,]  8.969159e+00
[10,]  8.840289e+00
[11,] -2.999143e+00
[12,] -2.777280e+00
[13,] -5.970016e+00
[14,] -6.063008e+00
[15,]  1.175424e+01
[16,]  1.165759e+01
[17,] -1.625296e+01
[18,] -1.582351e+01
[19,] -1.731430e+01
[20,] -1.571913e+01
[21,]  8.690474e+00
[22,]  8.137645e+00
[23,]  2.270214e+01
[24,]  2.093069e+01
> 
> # levels of fixed effect 1 in traits 1
> lvls_t1_f1 = rep(0, lvls_hys1)
> for ( i in 1 : lvls_hys1){
+     if (i == 1){
+         lvls_t1_f1[i] = 1
+     } else {
+         lvls_t1_f1[i] = 1 + (i - 1) * no_trts
+     }
+ }
> 
> # print
> lvls_t1_f1
[1] 1 3
> 
> # levels of fixed effect 1 in traits 2
> lvls_t2_f1 = lvls_t1_f1 + 1
> 
> # print
> lvls_t2_f1
[1] 2 4
> 
> # levels of fixed effect 2 in traits 1
> lvls_t1_f2 = rep(0, lvls_hys2)
> for ( i in 1 : lvls_hys2){
+     if (i == 1){
+         lvls_t1_f2[i] = 1
+     } else {
+         lvls_t1_f2[i] = 1 + (i - 1) * no_trts
+     }
+ }
> 
> # print
> lvls_t1_f2 = lvls_t1_f2 + no_trts * lvls_hys1
> lvls_t1_f2
[1] 5 7
> 
> # levels of fixed effect 2 in traits 2
> lvls_t2_f2 = lvls_t1_f2 + 1
> 
> # print
> lvls_t2_f2
[1] 6 8
> 
> # levels of animal effect in traits 1
> lvls_t1_ani = rep(0, lvls_ani)
> for ( i in 1 : lvls_ani){
+     if (i == 1){
+         lvls_t1_ani[i] = 1
+     } else {
+         lvls_t1_ani[i] = 1 + (i - 1) * no_trts
+     }
+ }
> lvls_t1_ani = lvls_t1_ani + no_trts * (lvls_hys1 + lvls_hys2)
> 
> # print
> lvls_t1_ani
[1]  9 11 13 15 17 19 21 23
> 
> # levels of fixed effect in traits 1
> lvls_t2_ani = lvls_t1_ani + 1
> 
> # print
> lvls_t2_ani
[1] 10 12 14 16 18 20 22 24
> 
> # solutions for fixed effects 1
> sol_t1_f1 = as.matrix(sol[lvls_t1_f1])
> sol_t2_f1 = as.matrix(sol[lvls_t2_f1])
> sol_f1 = cbind(sol_t1_f1, sol_t2_f1)
> 
> # print
> sol_f1
         [,1]         [,2]
[1,] 175.7313 1.532129e-12
[2,] 219.6133 1.427759e-12
> 
> # solutions for fixed effects 2
> sol_t1_f2 = as.matrix(sol[lvls_t1_f2])
> sol_t2_f2 = as.matrix(sol[lvls_t2_f2])
> sol_f2 = cbind(sol_t1_f2, sol_t2_f2)
> 
> # print
> sol_f2
              [,1]     [,2]
[1,] -6.947473e-14 243.2391
[2,]  6.444061e-13 240.5497
> 
> # breedinv value
> sol_t1_bv = sol[lvls_t1_ani]
> sol_t2_bv = sol[lvls_t2_ani]
> sol_bv = cbind(sol_t1_bv, sol_t2_bv)
> 
> # print
> sol_bv
      sol_t1_bv  sol_t2_bv
[1,]   8.969159   8.840289
[2,]  -2.999143  -2.777280
[3,]  -5.970016  -6.063008
[4,]  11.754242  11.657588
[5,] -16.252957 -15.823508
[6,] -17.314297 -15.719126
[7,]   8.690474   8.137645
[8,]  22.702139  20.930688
> 
> # reliability(r2), accuracy(r), standard error of prediction(SEP)
> 
> # diagonal elements of the generalized inverse of LHS for animal equation
> d_ani_t1 = diag(gi_lhs[lvls_t1_ani, lvls_t1_ani])
> d_ani_t2 = diag(gi_lhs[lvls_t2_ani, lvls_t2_ani])
> d_ani = cbind(d_ani_t1, d_ani_t2)
> 
> # print
> d_ani
     d_ani_t1 d_ani_t2
[1,] 32.09312 27.31163
[2,] 33.68775 28.74419
[3,] 32.33176 27.84285
[4,] 30.44248 26.38380
[5,] 29.58983 25.58931
[6,] 31.16542 26.56940
[7,] 30.77347 26.50225
[8,] 33.45367 28.61921
> 
> # reliability
> rel = matrix(rep(0, lvls_ani * no_trts), ncol = no_trts)
> 
> for (i in 1 : lvls_ani) {
+     rel[i, ] = 1 - d_ani[i, ]/diag(G)
+ }
> 
> # print
> rel
           [,1]       [,2]
[1,] 0.08305384 0.08961238
[2,] 0.03749289 0.04186017
[3,] 0.07623551 0.07190492
[4,] 0.13021474 0.12053983
[5,] 0.15457617 0.14702313
[6,] 0.10955930 0.11435322
[7,] 0.12075790 0.11659162
[8,] 0.04418092 0.04602625
> 
> # accuracy
> acc = sqrt(rel)
> 
> # print
> acc
          [,1]      [,2]
[1,] 0.2881906 0.2993533
[2,] 0.1936308 0.2045976
[3,] 0.2761078 0.2681509
[4,] 0.3608528 0.3471885
[5,] 0.3931618 0.3834359
[6,] 0.3309974 0.3381615
[7,] 0.3475024 0.3414551
[8,] 0.2101926 0.2145373
> 
> # standard error of prediction(SEP)
> sep = sqrt(d_ani)
> 
> # 확인
> sep
     d_ani_t1 d_ani_t2
[1,] 5.665079 5.226053
[2,] 5.804115 5.361361
[3,] 5.686102 5.276633
[4,] 5.517471 5.136517
[5,] 5.439654 5.058587
[6,] 5.582600 5.154552
[7,] 5.547384 5.148034
[8,] 5.783915 5.349693
> 
> # partitioning of breeding values
> 
> # yield deviation
> yd1 = matrix(rep(0, lvls_ani * no_trts), ncol = no_trts)
> 
> # numerator of n2
> a2 = array(rep(0, no_trts * no_trts * lvls_ani), dim = c(no_trts, no_trts, lvls_ani))
> 
> # looping data
> for (i in 1:no_data) {
+     yd1[data[i, 1], ] = as.matrix(data2[i, c(4, 5)] - sol_f1[data2[i, 2], ] - sol_f2[data2[i, 3], ])
+     
+     a2[, , data[i, 1]] = Ri[, , data2$dtsts[i]]
+ }
> 
> # print
> yd1
          [,1]       [,2]
[1,]   0.00000   0.000000
[2,]   0.00000   0.000000
[3,]   0.00000   0.000000
[4,]  25.26873  36.760913
[5,] -25.73127 -40.549726
[6,] -59.61329 -53.239087
[7,]   4.26873   6.760913
[8,]  65.38671  59.450274
> a2
, , 1

     [,1] [,2]
[1,]    0    0
[2,]    0    0

, , 2

     [,1] [,2]
[1,]    0    0
[2,]    0    0

, , 3

     [,1] [,2]
[1,]    0    0
[2,]    0    0

, , 4

             [,1]         [,2]
[1,]  0.018319812 -0.007066213
[2,] -0.007066213  0.017011254

, , 5

             [,1]         [,2]
[1,]  0.018319812 -0.007066213
[2,] -0.007066213  0.017011254

, , 6

             [,1]         [,2]
[1,]  0.018319812 -0.007066213
[2,] -0.007066213  0.017011254

, , 7

             [,1]         [,2]
[1,]  0.018319812 -0.007066213
[2,] -0.007066213  0.017011254

, , 8

             [,1]         [,2]
[1,]  0.018319812 -0.007066213
[2,] -0.007066213  0.017011254

> 
> # Parents average, progeny contribution
> 
> # parents avearge
> pa1 = matrix(rep(0, lvls_ani * no_trts), ncol = no_trts)
> 
> # progeny contribution numerator
> pc0 = matrix(rep(0, lvls_ani * no_trts), ncol = no_trts)
> 
> # numerator of n3, denominator of progeny contribution
> a3 = array(rep(0, no_trts * no_trts * lvls_ani), dim = c(no_trts, no_trts, lvls_ani))
> 
> # numerator of n1
> a1 = array(rep(0, no_trts * no_trts * lvls_ani), dim = c(no_trts, no_trts, lvls_ani))
> 
> 
> # looping pedi
> for (i in 1 : lvls_ani) {
+     
+     sire = pedi[i, 2]
+     dam = pedi[i, 3]
+     
+     if (sire == 0 & dam == 0) {
+         # both parents unknown PA
+         a1[, , i] = 1 * Gi
+         
+     } else if (sire != 0 & dam == 0) {
+         # 개체의 부만 알고 있을 경우
+         
+         # PA
+         a1[, , i] = 4/3 * Gi
+         pa1[i, ] = sol_bv[sire, ]/2
+         
+         # PC for sire
+         a3[, , sire] = a3[, , sire] + 0.5 * Gi * (2/3)
+         pc0[sire, ] = pc0[sire, ] + (0.5 * Gi * (2/3)) %*% (2 * sol_bv[i, ])
+         
+     } else if (sire == 0 & dam != 0) {
+         # dam known
+         
+         # PA
+         a1[, , i] = 4/3 * Gi
+         pa1[i, ] = sol_bv[dam, ]/2
+         
+         # PC for dam
+         a3[, , dam] = a3[, , dam] + 0.5 * Gi * (2/3)
+         pc0[dam, ] = pc0[dam, ] + (0.5 * Gi * (2/3)) %*% (2 * sol_bv[i])
+         
+     } else {
+         # both parents known
+         
+         # PA
+         a1[, , i] = 2 * Gi
+         pa1[i, ] = (sol_bv[sire, ] + sol_bv[dam, ])/2
+         
+         # PC for sire
+         a3[, , sire] = a3[, , sire] + 0.5 * Gi
+         pc0[sire, ] = pc0[sire, ] + (0.5 * Gi) %*% (2 * sol_bv[i, ] - sol_bv[dam, ])
+         
+         # PC for dam
+         a3[, , dam] = a3[, , dam] + 0.5 * Gi
+         pc0[dam, ] = pc0[dam, ] + (0.5 * Gi) %*% (2 * sol_bv[i, ] - sol_bv[sire, ])
+         
+     }
+ }
> 
> # print
> a1
, , 1

           [,1]       [,2]
[1,]  0.1127820 -0.1052632
[2,] -0.1052632  0.1315789

, , 2

           [,1]       [,2]
[1,]  0.1127820 -0.1052632
[2,] -0.1052632  0.1315789

, , 3

           [,1]       [,2]
[1,]  0.1127820 -0.1052632
[2,] -0.1052632  0.1315789

, , 4

           [,1]       [,2]
[1,]  0.2255639 -0.2105263
[2,] -0.2105263  0.2631579

, , 5

           [,1]       [,2]
[1,]  0.2255639 -0.2105263
[2,] -0.2105263  0.2631579

, , 6

           [,1]       [,2]
[1,]  0.2255639 -0.2105263
[2,] -0.2105263  0.2631579

, , 7

           [,1]       [,2]
[1,]  0.2255639 -0.2105263
[2,] -0.2105263  0.2631579

, , 8

           [,1]       [,2]
[1,]  0.2255639 -0.2105263
[2,] -0.2105263  0.2631579

> pa1
          [,1]      [,2]
[1,]  0.000000  0.000000
[2,]  0.000000  0.000000
[3,]  0.000000  0.000000
[4,]  2.985008  3.031504
[5,] -4.484580 -4.420144
[6,] -3.641899 -3.491610
[7,]  2.892113  2.797290
[8,]  8.829816  8.488967
> a3
, , 1

           [,1]       [,2]
[1,]  0.1691729 -0.1578947
[2,] -0.1578947  0.1973684

, , 2

           [,1]       [,2]
[1,]  0.1127820 -0.1052632
[2,] -0.1052632  0.1315789

, , 3

           [,1]       [,2]
[1,]  0.1127820 -0.1052632
[2,] -0.1052632  0.1315789

, , 4

            [,1]        [,2]
[1,]  0.05639098 -0.05263158
[2,] -0.05263158  0.06578947

, , 5

            [,1]        [,2]
[1,]  0.05639098 -0.05263158
[2,] -0.05263158  0.06578947

, , 6

     [,1] [,2]
[1,]    0    0
[2,]    0    0

, , 7

            [,1]        [,2]
[1,]  0.05639098 -0.05263158
[2,] -0.05263158  0.06578947

, , 8

     [,1] [,2]
[1,]    0    0
[2,]    0    0

> pc0
            [,1]        [,2]
[1,]  0.20250651  0.54768464
[2,] -0.09180779 -0.09946475
[3,] -0.07019742 -0.33868297
[4,]  0.14108377  0.24062679
[5,] -0.33859671 -0.35528541
[6,]  0.00000000  0.00000000
[7,]  0.31666002  0.25480212
[8,]  0.00000000  0.00000000
> 
> # denominator of n1, n2, n3, diagonal of animals in LHS
> n_de = a1 + a2 + a3
> 
> # print
> n_de
, , 1

           [,1]       [,2]
[1,]  0.2819549 -0.2631579
[2,] -0.2631579  0.3289474

, , 2

           [,1]       [,2]
[1,]  0.2255639 -0.2105263
[2,] -0.2105263  0.2631579

, , 3

           [,1]       [,2]
[1,]  0.2255639 -0.2105263
[2,] -0.2105263  0.2631579

, , 4

           [,1]       [,2]
[1,]  0.3002747 -0.2702241
[2,] -0.2702241  0.3459586

, , 5

           [,1]       [,2]
[1,]  0.3002747 -0.2702241
[2,] -0.2702241  0.3459586

, , 6

           [,1]       [,2]
[1,]  0.2438837 -0.2175925
[2,] -0.2175925  0.2801691

, , 7

           [,1]       [,2]
[1,]  0.3002747 -0.2702241
[2,] -0.2702241  0.3459586

, , 8

           [,1]       [,2]
[1,]  0.2438837 -0.2175925
[2,] -0.2175925  0.2801691

> 
> # parents average fraction of breeding values
> pa = matrix(rep(0, lvls_ani * no_trts), ncol = no_trts)
> for (i in 1 : lvls_ani) {
+     pa[i, ] = ginv(n_de[, , i]) %*% a1[, , i] %*% pa1[i, ]
+ }
> 
> # print
> pa
          [,1]      [,2]
[1,]  0.000000  0.000000
[2,]  0.000000  0.000000
[3,]  0.000000  0.000000
[4,]  1.876213  1.954972
[5,] -2.826254 -2.840790
[6,] -2.731376 -2.664305
[7,]  1.825628  1.793832
[8,]  6.620314  6.480231
> 
> # yield deviation fraction of breeding values
> yd = matrix(rep(0, lvls_ani * no_trts), ncol = no_trts)
> for (i in 1 : lvls_ani) {
+     yd[i, ] = ginv(n_de[, , i]) %*% a2[, , i] %*% yd1[i, ]
+ }
> 
> # print
> yd
           [,1]       [,2]
[1,]   0.000000   0.000000
[2,]   0.000000   0.000000
[3,]   0.000000   0.000000
[4,]   6.189546   6.126050
[5,]  -6.520144  -6.561129
[6,] -14.582921 -13.054821
[7,]   1.084030   1.091977
[8,]  16.081825  14.450457
> 
> # progeny contribution
> pc1 = matrix(rep(0, lvls_ani * no_trts), ncol = no_trts)
> for (i in 1 : lvls_ani) {
+     pc1[i, ] = ginv(a3[, , i]) %*% pc0[i, ]
+ }
> pc1[is.nan(pc1) == TRUE] = 0
> pc1
           [,1]       [,2]
[1,]  14.948599  14.733814
[2,]  -5.998286  -5.554561
[3,] -11.940033 -12.126017
[4,]  23.350964  22.338298
[5,] -43.597753 -40.278541
[6,]   0.000000   0.000000
[7,]  36.435120  33.021088
[8,]   0.000000   0.000000
> 
> # progeny contribution fraction of breeding values
> pc = matrix(rep(0, lvls_ani * no_trts), ncol = no_trts)
> for (i in 1 : lvls_ani) {
+     pc[i, ] = ginv(n_de[, , i]) %*% a3[, , i] %*% pc1[i, ]
+ }
> 
> # print
> pc
          [,1]      [,2]
[1,]  8.969159  8.840289
[2,] -2.999143 -2.777280
[3,] -5.970016 -6.063008
[4,]  3.688483  3.576566
[5,] -6.906559 -6.421589
[6,]  0.000000  0.000000
[7,]  5.780815  5.251836
[8,]  0.000000  0.000000
> 
> # Progeny(Daughter) Yield Deviation(PYD, DYD)
> 
> # n2 of progeny
> n2prog = array(rep(0, no_trts * no_trts * lvls_ani), dim = c(no_trts, no_trts, lvls_ani))
> for (i in 1 : lvls_ani) {
+     n2prog[, , i] = ginv((a1 + a2)[, , i]) %*% a2[, , i]
+ }
> 
> # print
> n2prog
, , 1

     [,1] [,2]
[1,]    0    0
[2,]    0    0

, , 2

     [,1] [,2]
[1,]    0    0
[2,]    0    0

, , 3

     [,1] [,2]
[1,]    0    0
[2,]    0    0

, , 4

          [,1]       [,2]
[1,] 0.1713401 0.08205946
[2,] 0.1078496 0.12444903

, , 5

          [,1]       [,2]
[1,] 0.1713401 0.08205946
[2,] 0.1078496 0.12444903

, , 6

          [,1]       [,2]
[1,] 0.1713401 0.08205946
[2,] 0.1078496 0.12444903

, , 7

          [,1]       [,2]
[1,] 0.1713401 0.08205946
[2,] 0.1078496 0.12444903

, , 8

          [,1]       [,2]
[1,] 0.1713401 0.08205946
[2,] 0.1078496 0.12444903

> 
> # numerator of dyd : summation of u of progeny * n2 of progeny * (2 * YD - bv of mate)
> dyd_n = matrix(rep(0, lvls_ani * no_trts), ncol = no_trts)
> # denominator of dyd : summation of u of progeny * n2 of progeny
> dyd_d = array(rep(0, no_trts * no_trts * lvls_ani), dim = c(no_trts, no_trts, lvls_ani))
> 
> # looping pedi
> for (i in 1 : lvls_ani) {
+     # i = 5
+     sire = pedi[i, 2]
+     dam = pedi[i, 3]
+     
+     if (sire == 0 & dam == 0) {
+         # both parents unknown
+         
+     } else if (sire != 0 & dam == 0) {
+         # 개체의 부만 알고 있을 경우
+         
+         # dyd_n
+         dyd_n[sire, ] = dyd_n[sire, ] + (n2prog[, , i] * 2/3) %*% (2 * yd1[i, ])
+         # dyd_d
+         dyd_d[, , sire] = dyd_d[, , sire] + n2prog[, , i] * 2/3
+         
+     } else if (sire == 0 & dam != 0) {
+         # dam known
+         
+         # dyd_n
+         dyd_n[dam, ] = dyd_n[dam, ] + (n2prog[, , i] * 2/3) %*% (2 * yd1[i, ])
+         # dyd_d
+         dyd_d[, , dam] = dyd_d[, , dam] + n2prog[, , i] * 2/3
+         
+     } else {
+         # both parents known
+         
+         # dyd_n
+         dyd_n[sire, ] = dyd_n[sire, ] + n2prog[, , i] %*% (2 * yd1[i, ] - sol_bv[dam, ])
+         dyd_n[dam, ] = dyd_n[dam, ] + n2prog[, , i] %*% (2 * yd1[i, ] - sol_bv[sire, ])
+         
+         # dyd_d
+         dyd_d[, , sire] = dyd_d[, , sire] + n2prog[, , i]
+         dyd_d[, , dam] = dyd_d[, , dam] + n2prog[, , i]
+         
+     }
+ }
> 
> # print
> dyd_n
           [,1]       [,2]
[1,]  20.358295  19.832631
[2,]  -1.522095  -1.711884
[3,] -15.128986 -15.088801
[4,]   4.092834   4.001939
[5,] -31.428049 -28.177127
[6,]   0.000000   0.000000
[7,]  29.901444  26.833429
[8,]   0.000000   0.000000
> dyd_d
, , 1

          [,1]      [,2]
[1,] 0.5140204 0.2461784
[2,] 0.3235487 0.3733471

, , 2

          [,1]      [,2]
[1,] 0.3426803 0.1641189
[2,] 0.2156991 0.2488981

, , 3

          [,1]      [,2]
[1,] 0.3426803 0.1641189
[2,] 0.2156991 0.2488981

, , 4

          [,1]       [,2]
[1,] 0.1713401 0.08205946
[2,] 0.1078496 0.12444903

, , 5

          [,1]       [,2]
[1,] 0.1713401 0.08205946
[2,] 0.1078496 0.12444903

, , 6

     [,1] [,2]
[1,]    0    0
[2,]    0    0

, , 7

          [,1]       [,2]
[1,] 0.1713401 0.08205946
[2,] 0.1078496 0.12444903

, , 8

     [,1] [,2]
[1,]    0    0
[2,]    0    0

> 
> # dyd
> dyd = matrix(rep(0, lvls_ani * no_trts), ncol = no_trts)
> 
> # looping pedi
> for (i in 1 : lvls_ani) {
+     dyd[i, ] = ginv(dyd_d[, , i]) %*% dyd_n[i, ]
+ }
> dyd[is.nan(dyd) == TRUE] = 0
> 
> # print
> dyd
            [,1]        [,2]
[1,]   24.215303   32.135781
[2,]   -1.962111   -5.177453
[3,]  -25.840090  -38.228967
[4,]   14.507476   19.584835
[5,] -128.195747 -115.318462
[6,]    0.000000    0.000000
[7,]  121.804253  110.060258
[8,]    0.000000    0.000000
> 
> # breeding values and fractions
> mt3_result = data.frame(animal = pedi[, 1], animal_bv = sol_bv, rel = rel, acc = acc, sep = sep, pa = pa, 
+     yd = yd, pc = pc, sum_of_fr = pa + yd + pc, dyd = dyd)
> 
> # print
> mt3_result
  animal animal_bv.sol_t1_bv animal_bv.sol_t2_bv      rel.1      rel.2     acc.1     acc.2 sep.d_ani_t1 sep.d_ani_t2
1      1            8.969159            8.840289 0.08305384 0.08961238 0.2881906 0.2993533     5.665079     5.226053
2      2           -2.999143           -2.777280 0.03749289 0.04186017 0.1936308 0.2045976     5.804115     5.361361
3      3           -5.970016           -6.063008 0.07623551 0.07190492 0.2761078 0.2681509     5.686102     5.276633
4      4           11.754242           11.657588 0.13021474 0.12053983 0.3608528 0.3471885     5.517471     5.136517
5      5          -16.252957          -15.823508 0.15457617 0.14702313 0.3931618 0.3834359     5.439654     5.058587
6      6          -17.314297          -15.719126 0.10955930 0.11435322 0.3309974 0.3381615     5.582600     5.154552
7      7            8.690474            8.137645 0.12075790 0.11659162 0.3475024 0.3414551     5.547384     5.148034
8      8           22.702139           20.930688 0.04418092 0.04602625 0.2101926 0.2145373     5.783915     5.349693
       pa.1      pa.2       yd.1       yd.2      pc.1      pc.2 sum_of_fr.1 sum_of_fr.2       dyd.1       dyd.2
1  0.000000  0.000000   0.000000   0.000000  8.969159  8.840289    8.969159    8.840289   24.215303   32.135781
2  0.000000  0.000000   0.000000   0.000000 -2.999143 -2.777280   -2.999143   -2.777280   -1.962111   -5.177453
3  0.000000  0.000000   0.000000   0.000000 -5.970016 -6.063008   -5.970016   -6.063008  -25.840090  -38.228967
4  1.876213  1.954972   6.189546   6.126050  3.688483  3.576566   11.754242   11.657588   14.507476   19.584835
5 -2.826254 -2.840790  -6.520144  -6.561129 -6.906559 -6.421589  -16.252957  -15.823508 -128.195747 -115.318462
6 -2.731376 -2.664305 -14.582921 -13.054821  0.000000  0.000000  -17.314297  -15.719126    0.000000    0.000000
7  1.825628  1.793832   1.084030   1.091977  5.780815  5.251836    8.690474    8.137645  121.804253  110.060258
8  6.620314  6.480231  16.081825  14.450457  0.000000  0.000000   22.702139   20.930688    0.000000    0.000000
> 
> # 파일로 출력, 분리자는 ',', 따옴표 없고
> output_filename = gsub("[ -]", "", paste("mt3_result_", Sys.Date(), ".csv"))
> write.table(mt3_result, output_filename, sep = ",", row.names = FALSE, quote = FALSE)
> 

 

출력한 결과는 다음과 같다.

animal,animal_bv.sol_t1_bv,animal_bv.sol_t2_bv,rel.1,rel.2,acc.1,acc.2,sep.d_ani_t1,sep.d_ani_t2,pa.1,pa.2,yd.1,yd.2,pc.1,pc.2,sum_of_fr.1,sum_of_fr.2,dyd.1,dyd.2
1,8.96915914423549,8.84028862908404,0.0830538368044876,0.0896123767117357,0.288190625809528,0.299353264073963,5.66507861479812,5.22605287943472,0,0,0,0,8.96915914423558,8.84028862908413,8.96915914423558,8.84028862908413,24.2153032651908,32.1357811454664
2,-2.99914278847626,-2.77728027471879,0.0374928900538841,0.0418601748683809,0.193630808638202,0.204597592528311,5.80411482037649,5.36136127806629,0,0,0,0,-2.99914278847663,-2.77728027471909,-2.99914278847663,-2.77728027471909,-1.96211132151522,-5.17745334586381
3,-5.97001635575898,-6.0630083543651,0.0762355143663486,0.0719049213793291,0.276107794830839,0.268150930222756,5.68610209169496,5.27663267232239,0,0,0,0,-5.97001635575903,-6.06300835436521,-5.97001635575903,-6.06300835436521,-25.8400897487149,-38.2289668542267
4,11.7542424313521,11.6575875661634,0.130214739810671,0.120539834590674,0.360852795209724,0.34718847128134,5.51747080704796,5.13651681222595,1.87621319870556,1.95497229270165,6.18954647427461,6.12604955346034,3.68848275837195,3.5765657200015,11.7542424313521,11.6575875661635,14.5074764284827,19.5848348700456
5,-16.2529566140671,-15.8235079782421,0.154576169145152,0.147023126540452,0.393161759515282,0.383435948419619,5.43965385662725,5.05858737235865,-2.82625401312119,-2.84079009320618,-6.52014354944942,-6.56112862579754,-6.90655905149658,-6.42158925923857,-16.2529566140672,-15.8235079782423,-128.195747122089,-115.318462113403
6,-17.3142968933337,-15.7191260030786,0.109559300197588,0.114353217167329,0.330997432312681,0.338161525261714,5.58260015522198,5.15455172493013,-2.7313756688042,-2.66430501516617,-14.5829212245295,-13.0548209879125,0,0,-17.3142968933337,-15.7191260030787,0,0
7,8.6904737239869,8.13764491523405,0.120757895768608,0.116591615294264,0.347502368004317,0.341455143897796,5.54738439700177,5.14803375485943,1.82562796263739,1.79383194633122,1.08403037069524,1.09197739220936,5.78081539065414,5.25183557669329,8.69047372398677,8.13764491523387,121.80425287791,110.060258438227
8,22.7021394832924,20.930688340762,0.0441809184024876,0.0460262510910151,0.210192574565534,0.214537295338165,5.78391457889144,5.34969274512748,6.62031436614155,6.48023116249169,16.0818251171507,14.4504571782701,0,0,22.7021394832923,20.9306883407618,0,0

 

csv 파일 이므로 엑셀에서 열면 다음과 같다.

+ Recent posts