Let's Study Animal Breeding

전체 글

pgAdmin 4를 이용한 Backup과 Restore

투정이 2020. 7. 27. 22:26

2020. 7. 27. 22:26

Restore를 하기 위하여 열심히 Backup을 하는데 정작 Restore는 그리 자주 할 일이 없어 자꾸 까먹는다. 자주 하면 텍스트 명령어로 하면 좋겠지만 자주 안하는데 굳이 텍스트 명령어로 할 일이 없을 것 같아 pgAdmin 4를 이용한 Backup과 Restore를 정리하고 필요할 때 보자.

중간에도 설명하지만 Backup은 텍스트 명령어로 매일 해야 한다.

pgadmin 실행

소유자가 myuser2 인 데이터베이스 mydb2

mydb2를 마우스 오른쪽 버튼 클릭 -> Backup 선택

백업 파일을 저장할 폴더를 지정하고, 파일 이름을 입력하기 위하여 ... 클릭

temp 폴더에 mydb2_backup 이란 이름으로 저장. 백업 포맷을 backup로 지정. create 클릭

format 은 custom, encoding은 utf-8 선택

Dump options 선택

위와 같이 설정하고 Backup을 클릭한다.

성공했다는 메시지가 나온다. More details ... 클릭하여 자세한 사항을 살펴 본다.

실제로 실행된 명령어와 백업 과정을 살펴볼 수 있다.

실제로 실행된 명령어는 다음과 같다.

C:\Program Files\PostgreSQL\12\bin\pg_dump.exe —file "C:\\Temp\\mydb2_backup.backup" --host "localhost" --port "5432" --username "postgres" --no-password —verbose --format=c --section=pre-data --section=data --section=post-data --column-inserts "mydb2"

백업은 자주 해야 하므로 위 명령을 배치 파일로 만들고 배치 파일을 작업 스케줄러로 주기적으로 실행시킬 수 있을 것이다.

PS. 2021.4.17.) 최근에 컴퓨터를 교체하면서 위의 백업 명령어를 cmd에서 실행하였더니, "no password supplied" 때문에 접속을 못한다는 에러가 떴다. 분명히 pgAdmin4에서는 잘 실행되는 명령어가 cmd에서 안 되다니. 그건 pgAdmin4에서는 접속할 때 password를 입력하고 그걸 이용해서 pg_dump를 실행하는데 cmd에서는 password가 없어서 실행을 못한다는 뜻이다. 찾아보니 pgpass 파일이 필요하다. linux 에서는 .pgpass 파일이고 윈도우에서는 pgpass.conf 파일이다. pgpass.conf 파일이 있어야할 위치는 %appdata%postgresql 폴더이다. 윈도우즈 탐색기에서 %appdata% 라고 치면 적당한 폴더로 이동한다. postgresql 폴더가 없으면 postgresql 폴더를 만들고 pgpass.conf 파일을 만든다. 예를 들어 사용자가 abcd 이면 다음 폴더에 pgpass.conf 파일을 만든다.C:\Users\abcd\AppData\Roaming\postgresql\pgpass.conf

pgpass.conf 파일의 내용은 다음과 같다.

localhost : 5432 : mydb : myuser : mypassword

구체적으로 써 보면 다음과 같다.

localhost:5432:*:postgres:12345678로컬호스트에 포트 번호 5432에 모든 데이터베이스에 사용자 postgres로 비밀번호 12345678로 접속하란 뜻이다. 이렇게 하고 postgresql 서버를 재시작하고 하니 pg_dump가 잘 동작한다. postgresql을 반드시 재시작해야 하는 지는 잘 모르겠다.

(2021.4.17.)

이제 새로운 컴퓨터에 리스토어 한다고 가정하자.

새로 유저와 빈 데이터베이스를 만들자.

새로운 Role과 Database가 만들어진 상태는 다음과 같다.

myuser2와 mydb2를 볼 수 있다. mydb2를 마우스 오른쪽 버튼으로 클릭하고 Restore ...를 클릭한다.

복원할 백업 파일을 선택하기 위하여 “...” 버튼을 클릭한다.

백업 파일을 선택하고, Format으로 backup을 선택하고 Select를 클릭한 후, Restore options를 클릭한다.

위와 같이 선택하고 Restore를 클릭한다.

성공적으로 완료되었음을 알 수 있다. More details ...를 클릭하여 자세히 알아 본다.

실제로 실행된 명령어를 확인할 수 있다.

다음은 리스토어된 데이터베이스와 테이블이다.

저작자표시

'ETC' 카테고리의 다른 글

레노보 아이디어패드 슬림3-15IIL05 구입 (2)	2020.10.27
레노보 아이디어패드 슬림3-17IML 3D 구입 (0)	2020.10.17
PostgreSQL psql을 이용한 role, database, table 생성 (0)	2020.07.27
The derivative of a determinant(행렬의 행렬식의 도함수, 미분) (1)	2020.05.21
Maximum Likelihood Estimation(MLE) (1)	2020.04.24

PostgreSQL psql을 이용한 role, database, table 생성

투정이 2020. 7. 27. 22:11

2020. 7. 27. 22:11

윈도우즈에서 postgresql을 사용하다 보니, pgAdmin과 EMS SQL Manager for PostgreSQL과 같은 그래픽 도구를 사용하므로 psql과 같은 텍스트 기반의 툴을 사용하는데 서툴었다.

간단히 정리하고 필요할 때마다 다시 보자.

psql을 다음과 같이 실행한다.

localhost에, postgres 데이터베이스에, 5432 포트로, postgres 사용자로 접속

###

bitnami wapp 로 설치할 경우 위와 같은 psql이 안 나올 수 있다. 그럴 경우 cmd에서 psql -U postgres 명령어로 psql에 접속

###

create user myuser2 with password '1234';

비밀번호 1234인 사용자 myuser2를 생성

create database mydb2 owner myuser2;

소유자가 myuser2인, mydb2 데이터베이스를 생성

select * from pg_database;

생성한 데이터베이스 확인.

\c mydb2 myuser2

mydb2에 접속

create table t1 (a varchar(2));

테이블 생성

select * from pg_tables where tableowner = 'myuser2';

생성한 테이블 확인

copy t1 from stdin;

표준 입력에서 테이블에 자료 입력

select * from t1;

입력한 자료 확인

select a, length(a) as len_a from t1 order by a;

a의 길이를 출력하고, a에 의해서 정렬

저작자표시

'ETC' 카테고리의 다른 글

레노보 아이디어패드 슬림3-17IML 3D 구입 (0)	2020.10.17
pgAdmin 4를 이용한 Backup과 Restore (0)	2020.07.27
The derivative of a determinant(행렬의 행렬식의 도함수, 미분) (1)	2020.05.21
Maximum Likelihood Estimation(MLE) (1)	2020.04.24
회귀분석 - Piecewise regression (0)	2020.03.01

R을 이용한 다형질 개체 모형(same model and missing record)의 육종가 구하기 프로그래밍

투정이 2020. 7. 26. 20:30

2020. 7. 26. 20:30

multiple trait animal model, same model and missing record

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

Raphael Mrode

Example 5.2 p78

자료

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

mt02_data.txt 로 저장

혈통

mt02_pedi.txt 로 저장

자료와 혈통을 읽어 육종가를 구하는 R 프로그램이다.

신뢰도, 정확도, Standard error of prediction을 구한다.

육종가를 PA(parenta average), YD(yield devation), PC(progeny contribution)으로 구한다.

개체의 DYD(daughter yield deviation)을 구한다.

# multiple trait animal model, same model and missing record - Date : 2020-07-26

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

# Raphael Mrode

# Example 5.2 p78

# 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/05_multiple_traits_02_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("mt02_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("mt02_data.txt", header = FALSE, sep = "", col.names = c("animal", "sex", "wwg", "pwg"), 
    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$wwg != 0, data2$dtsts + 1, data2$dtsts)
data2$dtsts = ifelse(data2$pwg != 0, data2$dtsts + 2, data2$dtsts)

# print
data2

# levels of animal and sex
lvls_ani = max(data2$animal)
lvls_sex = max(data2$sex)

# print
lvls_ani
lvls_sex

# variance component additive genetic
G = matrix(c(20, 18, 18, 40), 2, 2)
# residual
R = matrix(c(40, 11, 11, 30), 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_ani + lvls_sex))^2), nrow = no_trts * (lvls_ani + lvls_sex))

# print
dim(lhs)
lhs

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

# print
dim(rhs)
rhs

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

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

# print lhs and rhs
lhs
rhs

# inverse matrix of NRM
ai = ainv(pedi)

# print
ai

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

# print
lhs

# generalised inverse of lhs
gi_lhs = ginv(lhs)

# print
gi_lhs

# solution
sol = gi_lhs %*% rhs

# print
sol

# levels of fixed effect in traits 1
lvls_t1_f1 = rep(0, lvls_sex)
for ( i in 1 : lvls_sex){
    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 in traits 1
lvls_t2_f1 = lvls_t1_f1 + 1

# print
lvls_t2_f1

# 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_sex    

# 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 and animal effects
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

# 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(3, 4)] - sol_f1[data2[i, 2], ])
    
    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
mt2_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
mt2_result

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

다음은 실행 결과이다.

> # multiple trait animal model, same model and missing record - Date : 2020-07-26
> 
> # Linear Models for the Prediction of Animal Breeding Values, 3rd Edtion
> 
> # Raphael Mrode
> 
> # Example 5.2 p78
> 
> # 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 * lvls[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/05_multiple_traits_02_R")
> 
> # print working directory
> getwd()
[1] "D:/temp/05_multiple_traits_02_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("mt02_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   0
5      5    3   2
6      6    1   2
7      7    4   5
8      8    3   6
9      9    7   0
> 
> # read data : animal, dam, sex, weaning_weight
> data = read.table("mt02_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 0.0
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
6      9   2 4.0 0.0
> 
> # number of data
> no_data = nrow(data)
> 
> # print
> no_data
[1] 6
> 
> # how many traits does animal have?
> data2 = data.frame(data, dtsts = c(0))
> data2$dtsts = ifelse(data2$wwg != 0, data2$dtsts + 1, data2$dtsts)
> data2$dtsts = ifelse(data2$pwg != 0, data2$dtsts + 2, data2$dtsts)
> 
> # print
> data2
  animal sex wwg pwg dtsts
1      4   1 4.5 0.0     1
2      5   2 2.9 5.0     3
3      6   2 3.9 6.8     3
4      7   1 3.5 6.0     3
5      8   1 5.0 7.5     3
6      9   2 4.0 0.0     1
> 
> # levels of animal and sex
> lvls_ani = max(data2$animal)
> lvls_sex = max(data2$sex)
> 
> # print
> lvls_ani
[1] 9
> lvls_sex
[1] 2
> 
> # variance component additive genetic
> G = matrix(c(20, 18, 18, 40), 2, 2)
> # residual
> R = matrix(c(40, 11, 11, 30), 2, 2)
> 
> # print
> G
     [,1] [,2]
[1,]   20   18
[2,]   18   40
> R
     [,1] [,2]
[1,]   40   11
[2,]   11   30
> 
> # inverse of G
> Gi = ginv(G)
> 
> # print
> Gi
       [,1]   [,2]
[1,]  0.084 -0.038
[2,] -0.038  0.042
> 
> # 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.025    0
[2,] 0.000    0

, , 2

     [,1]  [,2]
[1,]    0 0.000
[2,]    0 0.033

, , 3

       [,1]   [,2]
[1,]  0.028 -0.010
[2,] -0.010  0.037

> 
> # empty lhs
> lhs = matrix(rep(0, (no_trts * (lvls_ani + lvls_sex))^2), nrow = no_trts * (lvls_ani + lvls_sex))
> 
> # print
> dim(lhs)
[1] 22 22
> lhs
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
 [1,]    0    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     0
 [3,]    0    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     0
 [5,]    0    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     0
 [7,]    0    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     0
 [9,]    0    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     0
[11,]    0    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     0
[13,]    0    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     0
[15,]    0    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     0
[17,]    0    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     0
[19,]    0    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     0
[21,]    0    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     0
> 
> # empty rhs
> rhs = matrix(rep(0, (no_trts * (lvls_ani + lvls_sex))), nrow = no_trts * (lvls_ani + lvls_sex))
> 
> # print
> dim(rhs)
[1] 22  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
> 
> # fill the MME
> for (i in 1:no_data) {
+ #i = 1
+     pos = pos_mme(no_trts, c(lvls_sex, lvls_ani), c(data2$sex[i], data2$animal[i]))
+ 
+     for (j in 1:length(pos)) {
+         for (k in 1:length(pos)) {
+             rfirst = pos[j]
+             rlast = (pos[j] + no_trts - 1)
+             cfirst = pos[k]
+             clast = (pos[k] + no_trts - 1)
+             
+             lhs[rfirst : rlast, cfirst : clast] = lhs[rfirst : rlast, cfirst : clast] + Ri[, , data2$dtsts[i]]
+         }
+         rhs[rfirst : rlast] = rhs[rfirst : rlast] + Ri[, , data2$dtsts[i]] %*% as.numeric(data2[i, 3:4])
+     }
+ }
> 
> # print lhs and rhs
> lhs
        [,1]   [,2]   [,3]   [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]  [,13]  [,14]  [,15]  [,16]  [,17]  [,18]  [,19]  [,20] [,21] [,22]
 [1,]  0.081 -0.020  0.000  0.000    0    0    0    0    0     0 0.025     0  0.000  0.000  0.000  0.000  0.028 -0.010  0.028 -0.010 0.000     0
 [2,] -0.020  0.074  0.000  0.000    0    0    0    0    0     0 0.000     0  0.000  0.000  0.000  0.000 -0.010  0.037 -0.010  0.037 0.000     0
 [3,]  0.000  0.000  0.081 -0.020    0    0    0    0    0     0 0.000     0  0.028 -0.010  0.028 -0.010  0.000  0.000  0.000  0.000 0.025     0
 [4,]  0.000  0.000 -0.020  0.074    0    0    0    0    0     0 0.000     0 -0.010  0.037 -0.010  0.037  0.000  0.000  0.000  0.000 0.000     0
 [5,]  0.000  0.000  0.000  0.000    0    0    0    0    0     0 0.000     0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.000     0
 [6,]  0.000  0.000  0.000  0.000    0    0    0    0    0     0 0.000     0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.000     0
 [7,]  0.000  0.000  0.000  0.000    0    0    0    0    0     0 0.000     0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.000     0
 [8,]  0.000  0.000  0.000  0.000    0    0    0    0    0     0 0.000     0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.000     0
 [9,]  0.000  0.000  0.000  0.000    0    0    0    0    0     0 0.000     0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.000     0
[10,]  0.000  0.000  0.000  0.000    0    0    0    0    0     0 0.000     0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.000     0
[11,]  0.025  0.000  0.000  0.000    0    0    0    0    0     0 0.025     0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.000     0
[12,]  0.000  0.000  0.000  0.000    0    0    0    0    0     0 0.000     0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.000     0
[13,]  0.000  0.000  0.028 -0.010    0    0    0    0    0     0 0.000     0  0.028 -0.010  0.000  0.000  0.000  0.000  0.000  0.000 0.000     0
[14,]  0.000  0.000 -0.010  0.037    0    0    0    0    0     0 0.000     0 -0.010  0.037  0.000  0.000  0.000  0.000  0.000  0.000 0.000     0
[15,]  0.000  0.000  0.028 -0.010    0    0    0    0    0     0 0.000     0  0.000  0.000  0.028 -0.010  0.000  0.000  0.000  0.000 0.000     0
[16,]  0.000  0.000 -0.010  0.037    0    0    0    0    0     0 0.000     0  0.000  0.000 -0.010  0.037  0.000  0.000  0.000  0.000 0.000     0
[17,]  0.028 -0.010  0.000  0.000    0    0    0    0    0     0 0.000     0  0.000  0.000  0.000  0.000  0.028 -0.010  0.000  0.000 0.000     0
[18,] -0.010  0.037  0.000  0.000    0    0    0    0    0     0 0.000     0  0.000  0.000  0.000  0.000 -0.010  0.037  0.000  0.000 0.000     0
[19,]  0.028 -0.010  0.000  0.000    0    0    0    0    0     0 0.000     0  0.000  0.000  0.000  0.000  0.000  0.000  0.028 -0.010 0.000     0
[20,] -0.010  0.037  0.000  0.000    0    0    0    0    0     0 0.000     0  0.000  0.000  0.000  0.000  0.000  0.000 -0.010  0.037 0.000     0
[21,]  0.000  0.000  0.025  0.000    0    0    0    0    0     0 0.000     0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.025     0
[22,]  0.000  0.000  0.000  0.000    0    0    0    0    0     0 0.000     0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 0.000     0
> rhs
       [,1]
 [1,] 0.211
 [2,] 0.414
 [3,] 0.169
 [4,] 0.368
 [5,] 0.000
 [6,] 0.000
 [7,] 0.000
 [8,] 0.000
 [9,] 0.000
[10,] 0.000
[11,] 0.112
[12,] 0.000
[13,] 0.030
[14,] 0.156
[15,] 0.039
[16,] 0.212
[17,] 0.036
[18,] 0.187
[19,] 0.063
[20,] 0.227
[21,] 0.100
[22,] 0.000
> 
> # inverse matrix of NRM
> ai = ainv(pedi)
> 
> # print
> ai
       [,1] [,2] [,3]  [,4] [,5] [,6]  [,7] [,8]  [,9]
 [1,]  1.83  0.5  0.0 -0.67  0.0 -1.0  0.00    0  0.00
 [2,]  0.50  2.0  0.5  0.00 -1.0 -1.0  0.00    0  0.00
 [3,]  0.00  0.5  2.0  0.00 -1.0  0.5  0.00   -1  0.00
 [4,] -0.67  0.0  0.0  1.83  0.5  0.0 -1.00    0  0.00
 [5,]  0.00 -1.0 -1.0  0.50  2.5  0.0 -1.00    0  0.00
 [6,] -1.00 -1.0  0.5  0.00  0.0  2.5  0.00   -1  0.00
 [7,]  0.00  0.0  0.0 -1.00 -1.0  0.0  2.33    0 -0.67
 [8,]  0.00  0.0 -1.0  0.00  0.0 -1.0  0.00    2  0.00
 [9,]  0.00  0.0  0.0  0.00  0.0  0.0 -0.67    0  1.33
> 
> # add ai to lhs
> afirst = no_trts * lvls_sex + 1
> alast = no_trts * (lvls_sex + lvls_ani)
> lhs[afirst : alast, afirst : alast] = lhs[afirst : alast, afirst : alast] + ai %x% Gi
> 
> # print
> lhs
        [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]  [,10]  [,11]  [,12]  [,13]  [,14]  [,15]  [,16]  [,17]  [,18]  [,19]  [,20]  [,21]  [,22]
 [1,]  0.081 -0.020  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.025  0.000  0.000  0.000  0.000  0.000  0.028 -0.010  0.028 -0.010  0.000  0.000
 [2,] -0.020  0.074  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 -0.010  0.037 -0.010  0.037  0.000  0.000
 [3,]  0.000  0.000  0.081 -0.020  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.028 -0.010  0.028 -0.010  0.000  0.000  0.000  0.000  0.025  0.000
 [4,]  0.000  0.000 -0.020  0.074  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 -0.010  0.037 -0.010  0.037  0.000  0.000  0.000  0.000  0.000  0.000
 [5,]  0.000  0.000  0.000  0.000  0.154 -0.069  0.042 -0.019  0.000  0.000 -0.056  0.025  0.000  0.000 -0.084  0.038  0.000  0.000  0.000  0.000  0.000  0.000
 [6,]  0.000  0.000  0.000  0.000 -0.069  0.077 -0.019  0.021  0.000  0.000  0.025 -0.028  0.000  0.000  0.038 -0.042  0.000  0.000  0.000  0.000  0.000  0.000
 [7,]  0.000  0.000  0.000  0.000  0.042 -0.019  0.168 -0.076  0.042 -0.019  0.000  0.000 -0.084  0.038 -0.084  0.038  0.000  0.000  0.000  0.000  0.000  0.000
 [8,]  0.000  0.000  0.000  0.000 -0.019  0.021 -0.076  0.084 -0.019  0.021  0.000  0.000  0.038 -0.042  0.038 -0.042  0.000  0.000  0.000  0.000  0.000  0.000
 [9,]  0.000  0.000  0.000  0.000  0.000  0.000  0.042 -0.019  0.168 -0.076  0.000  0.000 -0.084  0.038  0.042 -0.019  0.000  0.000 -0.084  0.038  0.000  0.000
[10,]  0.000  0.000  0.000  0.000  0.000  0.000 -0.019  0.021 -0.076  0.084  0.000  0.000  0.038 -0.042 -0.019  0.021  0.000  0.000  0.038 -0.042  0.000  0.000
[11,]  0.025  0.000  0.000  0.000 -0.056  0.025  0.000  0.000  0.000  0.000  0.179 -0.069  0.042 -0.019  0.000  0.000 -0.084  0.038  0.000  0.000  0.000  0.000
[12,]  0.000  0.000  0.000  0.000  0.025 -0.028  0.000  0.000  0.000  0.000 -0.069  0.077 -0.019  0.021  0.000  0.000  0.038 -0.042  0.000  0.000  0.000  0.000
[13,]  0.000  0.000  0.028 -0.010  0.000  0.000 -0.084  0.038 -0.084  0.038  0.042 -0.019  0.238 -0.105  0.000  0.000 -0.084  0.038  0.000  0.000  0.000  0.000
[14,]  0.000  0.000 -0.010  0.037  0.000  0.000  0.038 -0.042  0.038 -0.042 -0.019  0.021 -0.105  0.142  0.000  0.000  0.038 -0.042  0.000  0.000  0.000  0.000
[15,]  0.000  0.000  0.028 -0.010 -0.084  0.038 -0.084  0.038  0.042 -0.019  0.000  0.000  0.000  0.000  0.238 -0.105  0.000  0.000 -0.084  0.038  0.000  0.000
[16,]  0.000  0.000 -0.010  0.037  0.038 -0.042  0.038 -0.042 -0.019  0.021  0.000  0.000  0.000  0.000 -0.105  0.142  0.000  0.000  0.038 -0.042  0.000  0.000
[17,]  0.028 -0.010  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 -0.084  0.038 -0.084  0.038  0.000  0.000  0.224 -0.098  0.000  0.000 -0.056  0.025
[18,] -0.010  0.037  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.038 -0.042  0.038 -0.042  0.000  0.000 -0.098  0.135  0.000  0.000  0.025 -0.028
[19,]  0.028 -0.010  0.000  0.000  0.000  0.000  0.000  0.000 -0.084  0.038  0.000  0.000  0.000  0.000 -0.084  0.038  0.000  0.000  0.196 -0.086  0.000  0.000
[20,] -0.010  0.037  0.000  0.000  0.000  0.000  0.000  0.000  0.038 -0.042  0.000  0.000  0.000  0.000  0.038 -0.042  0.000  0.000 -0.086  0.121  0.000  0.000
[21,]  0.000  0.000  0.025  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 -0.056  0.025  0.000  0.000  0.137 -0.050
[22,]  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.025 -0.028  0.000  0.000 -0.050  0.056
> 
> # generalised inverse of lhs
> gi_lhs = ginv(lhs)
> 
> # print
> gi_lhs
       [,1]  [,2] [,3]  [,4]   [,5]    [,6]   [,7]    [,8]  [,9]  [,10]  [,11]  [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22]
 [1,]  23.8  13.1  6.1   5.6 -6.518  -5.796 -3.285  -2.956 -5.26  -4.88 -10.54  -9.34  -5.9  -5.4 -6.56  -5.9 -11.3 -10.0  -9.6  -8.8  -5.8  -5.1
 [2,]  13.1  38.4  6.1  15.4 -5.248 -10.861 -3.700  -9.280 -5.64 -13.84  -7.53 -14.27  -6.4 -15.9 -6.46 -15.0 -10.7 -24.4 -10.0 -24.0  -5.6 -12.4
 [3,]   6.1   6.1 22.5  12.0 -4.457  -4.020 -7.196  -6.477 -3.84  -3.50  -3.93  -3.45  -9.4  -8.4 -9.15  -8.3  -7.8  -6.9  -6.5  -6.0  -9.0  -8.0
 [4,]   5.6  15.4 12.0  37.6 -4.397  -9.923 -7.567 -18.710 -3.69  -9.18  -3.39  -6.15  -9.6 -23.5 -9.58 -23.4  -7.0 -15.4  -6.4 -15.9  -5.8  -9.8
 [5,]  -6.5  -5.2 -4.5  -4.4 18.567  16.171  0.049   0.044  1.47   1.90   8.71   7.50   1.9   2.5  8.26   6.7   5.6   5.2   5.2   4.7   3.3   3.0
 [6,]  -5.8 -10.9 -4.0  -9.9 16.171  36.405  0.024   0.021  1.90   3.78   7.55  17.54   2.5   5.1  6.63  15.1   5.2  11.2   4.7  10.0   3.0   6.0
 [7,]  -3.3  -3.7 -7.2  -7.6  0.049   0.024 19.055  17.149 -0.56  -0.45   0.91   0.76   8.8   8.0  9.01   8.1   4.9   4.3   4.1   3.8   3.7   3.3
 [8,]  -3.0  -9.3 -6.5 -18.7  0.044   0.021 17.149  39.234 -0.50  -0.40   0.82   0.68   8.0  19.1  8.11  19.2   4.4   9.8   3.7   9.3   3.4   6.0
 [9,]  -5.3  -5.6 -3.8  -3.7  1.466   1.901 -0.561  -0.505 17.93  15.39   2.14   2.46   7.6   6.0  0.98   1.6   5.2   4.9   8.5   7.3   2.9   2.7
[10,]  -4.9 -13.8 -3.5  -9.2  1.902   3.777 -0.450  -0.405 15.39  34.99   2.41   4.09   6.1  14.6  1.66   4.0   4.9  11.2   7.3  17.2   2.8   5.8
[11,] -10.5  -7.5 -3.9  -3.4  8.710   7.550  0.913   0.821  2.14   2.41  16.94  14.81   2.1   2.2  5.13   4.4   9.5   8.1   5.2   5.2   4.5   3.9
[12,]  -9.3 -14.3 -3.5  -6.1  7.505  17.544  0.759   0.683  2.46   4.09  14.81  35.43   2.2   3.3  4.26   8.9   8.0  17.4   5.2   9.7   3.9   8.5
[13,]  -5.9  -6.4 -9.4  -9.6  1.857   2.458  8.845   7.960  7.59   6.09   2.12   2.24  16.0  13.2  6.41   7.0   8.8   7.4   6.8   6.5   5.7   4.9
[14,]  -5.4 -15.9 -8.4 -23.5  2.459   5.077  7.974  19.076  6.05  14.58   2.24   3.30  13.2  31.3  7.00  16.6   7.5  16.9   6.4  15.7   5.0   9.6
[15,]  -6.6  -6.5 -9.2  -9.6  8.258   6.625  9.009   8.108  0.98   1.66   5.13   4.26   6.4   7.0 16.22  13.4   6.4   6.3   8.2   7.0   4.8   4.6
[16,]  -5.9 -15.0 -8.3 -23.4  6.666  15.057  8.112  19.201  1.59   4.00   4.38   8.89   7.0  16.6 13.36  31.3   6.3  13.9   6.9  16.6   4.6   8.2
[17,] -11.3 -10.7 -7.8  -7.0  5.614   5.168  4.877   4.389  5.18   4.95   9.51   8.01   8.8   7.5  6.39   6.3  17.2  14.2   7.2   7.8   8.4   6.9
[18,] -10.0 -24.4 -6.9 -15.4  5.194  11.222  4.295   9.816  4.87  11.16   8.12  17.37   7.4  16.9  6.26  13.9  14.2  32.0   7.6  17.3   7.0  15.9
[19,]  -9.6 -10.0 -6.5  -6.4  5.229   4.669  4.065   3.659  8.48   7.30   5.18   5.20   6.8   6.4  8.17   6.9   7.2   7.6  16.3  13.5   4.4   4.5
[20,]  -8.8 -24.0 -6.0 -15.9  4.700  10.018  3.758   9.332  7.28  17.20   5.17   9.67   6.5  15.7  7.05  16.6   7.8  17.3  13.5  31.7   4.5   9.1
[21,]  -5.8  -5.6 -9.0  -5.8  3.257   2.976  3.736   3.362  2.93   2.75   4.53   3.85   5.7   5.0  4.82   4.6   8.4   7.0   4.4   4.5  16.4  14.5
[22,]  -5.1 -12.4 -8.0  -9.8  3.002   5.964  3.316   5.959  2.74   5.83   3.86   8.55   4.9   9.6  4.59   8.2   6.9  15.9   4.5   9.1  14.5  35.8
> 
> # solution
> sol = gi_lhs %*% rhs
> 
> # print
> sol
        [,1]
 [1,]  4.367
 [2,]  6.834
 [3,]  3.657
 [4,]  6.007
 [5,]  0.130
 [6,]  0.266
 [7,] -0.084
 [8,] -0.075
 [9,] -0.098
[10,] -0.194
[11,]  0.007
[12,]  0.016
[13,] -0.343
[14,] -0.555
[15,]  0.192
[16,]  0.440
[17,] -0.308
[18,] -0.483
[19,]  0.201
[20,]  0.349
[21,] -0.018
[22,] -0.119
> 
> # levels of fixed effect in traits 1
> lvls_t1_f1 = rep(0, lvls_sex)
> for ( i in 1 : lvls_sex){
+     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 in traits 1
> lvls_t2_f1 = lvls_t1_f1 + 1
> 
> # print
> lvls_t2_f1
[1] 2 4
> 
> # 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_sex    
> 
> # print
> lvls_t1_ani
[1]  5  7  9 11 13 15 17 19 21
> 
> # levels of fixed effect in traits 1
> lvls_t2_ani = lvls_t1_ani + 1
> 
> # print
> lvls_t2_ani
[1]  6  8 10 12 14 16 18 20 22
> 
> # solutions for fixed effects and animal effects
> 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,]  4.4  6.8
[2,]  3.7  6.0
> 
> # 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,]     0.130     0.266
 [2,]    -0.084    -0.075
 [3,]    -0.098    -0.194
 [4,]     0.007     0.016
 [5,]    -0.343    -0.555
 [6,]     0.192     0.440
 [7,]    -0.308    -0.483
 [8,]     0.201     0.349
 [9,]    -0.018    -0.119
> 
> # 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,]       19       36
 [2,]       19       39
 [3,]       18       35
 [4,]       17       35
 [5,]       16       31
 [6,]       16       31
 [7,]       17       32
 [8,]       16       32
 [9,]       16       36
> 
> # 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.072 0.090
 [2,] 0.047 0.019
 [3,] 0.103 0.125
 [4,] 0.153 0.114
 [5,] 0.200 0.217
 [6,] 0.189 0.218
 [7,] 0.141 0.200
 [8,] 0.187 0.208
 [9,] 0.180 0.106
> 
> # accuracy
> acc = sqrt(rel)
> 
> # print
> acc
      [,1] [,2]
 [1,] 0.27 0.30
 [2,] 0.22 0.14
 [3,] 0.32 0.35
 [4,] 0.39 0.34
 [5,] 0.45 0.47
 [6,] 0.43 0.47
 [7,] 0.38 0.45
 [8,] 0.43 0.46
 [9,] 0.42 0.33
> 
> # standard error of prediction(SEP)
> sep = sqrt(d_ani)
> 
> # 확인
> sep
      d_ani_t1 d_ani_t2
 [1,]      4.3      6.0
 [2,]      4.4      6.3
 [3,]      4.2      5.9
 [4,]      4.1      6.0
 [5,]      4.0      5.6
 [6,]      4.0      5.6
 [7,]      4.1      5.7
 [8,]      4.0      5.6
 [9,]      4.1      6.0
> 
> # 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(3, 4)] - sol_f1[data2[i, 2], ])
+     
+     a2[, , data[i, 1]] = Ri[, , data2$dtsts[i]]
+ }
> 
> # print
> yd1
       [,1]  [,2]
 [1,]  0.00  0.00
 [2,]  0.00  0.00
 [3,]  0.00  0.00
 [4,]  0.13 -6.83
 [5,] -0.76 -1.01
 [6,]  0.24  0.79
 [7,] -0.87 -0.83
 [8,]  0.63  0.67
 [9,]  0.34 -6.01
> 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.025    0
[2,] 0.000    0

, , 5

       [,1]   [,2]
[1,]  0.028 -0.010
[2,] -0.010  0.037

, , 6

       [,1]   [,2]
[1,]  0.028 -0.010
[2,] -0.010  0.037

, , 7

       [,1]   [,2]
[1,]  0.028 -0.010
[2,] -0.010  0.037

, , 8

       [,1]   [,2]
[1,]  0.028 -0.010
[2,] -0.010  0.037

, , 9

      [,1] [,2]
[1,] 0.025    0
[2,] 0.000    0

> 
> # 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.084 -0.038
[2,] -0.038  0.042

, , 2

       [,1]   [,2]
[1,]  0.084 -0.038
[2,] -0.038  0.042

, , 3

       [,1]   [,2]
[1,]  0.084 -0.038
[2,] -0.038  0.042

, , 4

      [,1]   [,2]
[1,]  0.11 -0.050
[2,] -0.05  0.056

, , 5

       [,1]   [,2]
[1,]  0.168 -0.076
[2,] -0.076  0.084

, , 6

       [,1]   [,2]
[1,]  0.168 -0.076
[2,] -0.076  0.084

, , 7

       [,1]   [,2]
[1,]  0.168 -0.076
[2,] -0.076  0.084

, , 8

       [,1]   [,2]
[1,]  0.168 -0.076
[2,] -0.076  0.084

, , 9

      [,1]   [,2]
[1,]  0.11 -0.050
[2,] -0.05  0.056

> pa1
        [,1]   [,2]
 [1,]  0.000  0.000
 [2,]  0.000  0.000
 [3,]  0.000  0.000
 [4,]  0.065  0.133
 [5,] -0.091 -0.135
 [6,]  0.023  0.095
 [7,] -0.168 -0.269
 [8,]  0.047  0.123
 [9,] -0.154 -0.241
> a3
, , 1

       [,1]   [,2]
[1,]  0.070 -0.032
[2,] -0.032  0.035

, , 2

       [,1]   [,2]
[1,]  0.084 -0.038
[2,] -0.038  0.042

, , 3

       [,1]   [,2]
[1,]  0.084 -0.038
[2,] -0.038  0.042

, , 4

       [,1]   [,2]
[1,]  0.042 -0.019
[2,] -0.019  0.021

, , 5

       [,1]   [,2]
[1,]  0.042 -0.019
[2,] -0.019  0.021

, , 6

       [,1]   [,2]
[1,]  0.042 -0.019
[2,] -0.019  0.021

, , 7

       [,1]   [,2]
[1,]  0.028 -0.013
[2,] -0.013  0.014

, , 8

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

, , 9

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

> pc0
         [,1]     [,2]
 [1,]  0.0015  1.2e-02
 [2,] -0.0084 -4.1e-16
 [3,] -0.0018 -8.9e-03
 [4,] -0.0037 -3.5e-03
 [5,] -0.0076 -8.8e-03
 [6,]  0.0041  9.3e-03
 [7,]  0.0020 -2.9e-03
 [8,]  0.0000  0.0e+00
 [9,]  0.0000  0.0e+00
> 
> # denominator of n1, n2, n3, diagonal of animals in LHS
> n_de = a1 + a2 + a3
> 
> # print
> n_de
, , 1

       [,1]   [,2]
[1,]  0.154 -0.069
[2,] -0.069  0.077

, , 2

       [,1]   [,2]
[1,]  0.168 -0.076
[2,] -0.076  0.084

, , 3

       [,1]   [,2]
[1,]  0.168 -0.076
[2,] -0.076  0.084

, , 4

       [,1]   [,2]
[1,]  0.179 -0.069
[2,] -0.069  0.077

, , 5

      [,1]  [,2]
[1,]  0.24 -0.10
[2,] -0.10  0.14

, , 6

      [,1]  [,2]
[1,]  0.24 -0.10
[2,] -0.10  0.14

, , 7

       [,1]   [,2]
[1,]  0.224 -0.098
[2,] -0.098  0.135

, , 8

       [,1]   [,2]
[1,]  0.196 -0.086
[2,] -0.086  0.121

, , 9

      [,1]   [,2]
[1,]  0.14 -0.050
[2,] -0.05  0.056

> 
> # 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.000  0.000
 [2,]  0.000  0.000
 [3,]  0.000  0.000
 [4,]  0.037  0.088
 [5,] -0.052 -0.070
 [6,]  0.008  0.050
 [7,] -0.099 -0.146
 [8,]  0.025  0.074
 [9,] -0.112 -0.204
> 
> # 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.000  0.000
 [2,]  0.000  0.000
 [3,]  0.000  0.000
 [4,]  0.029  0.026
 [5,] -0.203 -0.358
 [6,]  0.115  0.274
 [7,] -0.208 -0.315
 [8,]  0.176  0.275
 [9,]  0.094  0.084
> 
> # 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,]  0.286  0.59
 [2,] -0.167 -0.15
 [3,] -0.196 -0.39
 [4,] -0.274 -0.41
 [5,] -0.623 -0.98
 [6,]  0.499  0.89
 [7,] -0.037 -0.24
 [8,]  0.000  0.00
 [9,]  0.000  0.00
> 
> # 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,]  0.12982  0.266
 [2,] -0.08361 -0.075
 [3,] -0.09807 -0.194
 [4,] -0.05861 -0.098
 [5,] -0.08802 -0.127
 [6,]  0.06827  0.116
 [7,] -0.00079 -0.022
 [8,]  0.00000  0.000
 [9,]  0.00000  0.000
> 
> # 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.27    0
[2,] 0.25    0

, , 5

      [,1] [,2]
[1,] 0.152 0.12
[2,] 0.024 0.39

, , 6

      [,1] [,2]
[1,] 0.152 0.12
[2,] 0.024 0.39

, , 7

      [,1] [,2]
[1,] 0.152 0.12
[2,] 0.024 0.39

, , 8

      [,1] [,2]
[1,] 0.152 0.12
[2,] 0.024 0.39

, , 9

     [,1] [,2]
[1,] 0.27    0
[2,] 0.25    0

> 
> # 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,]  0.33  0.71
 [2,] -0.22 -0.22
 [3,] -0.18 -0.42
 [4,] -0.34 -0.47
 [5,] -0.47 -0.70
 [6,]  0.39  0.63
 [7,]  0.12  0.11
 [8,]  0.00  0.00
 [9,]  0.00  0.00
> dyd_d
, , 1

     [,1] [,2]
[1,] 0.33 0.12
[2,] 0.19 0.39

, , 2

      [,1] [,2]
[1,] 0.305 0.24
[2,] 0.048 0.78

, , 3

      [,1] [,2]
[1,] 0.305 0.24
[2,] 0.048 0.78

, , 4

      [,1] [,2]
[1,] 0.152 0.12
[2,] 0.024 0.39

, , 5

      [,1] [,2]
[1,] 0.152 0.12
[2,] 0.024 0.39

, , 6

      [,1] [,2]
[1,] 0.152 0.12
[2,] 0.024 0.39

, , 7

     [,1] [,2]
[1,] 0.18    0
[2,] 0.16    0

, , 8

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

, , 9

     [,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,]  0.43  1.60
 [2,] -0.53 -0.25
 [3,] -0.18 -0.52
 [4,] -1.39 -1.11
 [5,] -1.74 -1.68
 [6,]  1.36  1.53
 [7,]  0.69  0.00
 [8,]  0.00  0.00
 [9,]  0.00  0.00
> 
> # breeding values and fractions
> mt2_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
> mt2_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   pa.1   pa.2   yd.1   yd.2     pc.1   pc.2 sum_of_fr.1
1      1               0.130               0.266 0.072 0.090  0.27  0.30          4.3          6.0  0.000  0.000  0.000  0.000  0.12982  0.266       0.130
2      2              -0.084              -0.075 0.047 0.019  0.22  0.14          4.4          6.3  0.000  0.000  0.000  0.000 -0.08361 -0.075      -0.084
3      3              -0.098              -0.194 0.103 0.125  0.32  0.35          4.2          5.9  0.000  0.000  0.000  0.000 -0.09807 -0.194      -0.098
4      4               0.007               0.016 0.153 0.114  0.39  0.34          4.1          6.0  0.037  0.088  0.029  0.026 -0.05861 -0.098       0.007
5      5              -0.343              -0.555 0.200 0.217  0.45  0.47          4.0          5.6 -0.052 -0.070 -0.203 -0.358 -0.08802 -0.127      -0.343
6      6               0.192               0.440 0.189 0.218  0.43  0.47          4.0          5.6  0.008  0.050  0.115  0.274  0.06827  0.116       0.192
7      7              -0.308              -0.483 0.141 0.200  0.38  0.45          4.1          5.7 -0.099 -0.146 -0.208 -0.315 -0.00079 -0.022      -0.308
8      8               0.201               0.349 0.187 0.208  0.43  0.46          4.0          5.6  0.025  0.074  0.176  0.275  0.00000  0.000       0.201
9      9              -0.018              -0.119 0.180 0.106  0.42  0.33          4.1          6.0 -0.112 -0.204  0.094  0.084  0.00000  0.000      -0.018
  sum_of_fr.2 dyd.1 dyd.2
1       0.266  0.43  1.60
2      -0.075 -0.53 -0.25
3      -0.194 -0.18 -0.52
4       0.016 -1.39 -1.11
5      -0.555 -1.74 -1.68
6       0.440  1.36  1.53
7      -0.483  0.69  0.00
8       0.349  0.00  0.00
9      -0.119  0.00  0.00
> 
> # 파일로 출력, 분리자는 ',', 따옴표 없고
> output_filename = gsub("[ -]", "", paste("mt2_result_", Sys.Date(), ".csv"))
> write.table(mt2_result, output_filename, sep = ",", row.names = FALSE, quote = FALSE)

출력 파일(mt2_result_20200723.csv)은 다음과 같다.

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,0.129823830715541,0.266245867684148,0.0716308723768238,0.0898698568225488,0.267639444732692,0.299783016234324,4.3089885765065,6.03367265660792,0,0,0,0,0.12982383071554,0.266245867684148,0.12982383071554,0.266245867684148,0.425630149209032,1.6034276646541
2,-0.0836108051833632,-0.0752497246650299,0.0472677904320411,0.0191434551249767,0.217411569223078,0.138359875415442,4.36516256184798,6.2637258716359,0,0,0,0,-0.0836108051833641,-0.0752497246650325,-0.0836108051833641,-0.0752497246650325,-0.529046133710049,-0.25000119324945
3,-0.0980748544120003,-0.194040392874268,0.103425883888995,0.125227278344733,0.321598948830674,0.353874664739838,4.23455810235497,5.91531139215939,0,0,0,0,-0.0980748544119998,-0.194040392874261,-0.0980748544119998,-0.194040392874261,-0.17737262403059,-0.523780147010833
4,0.0070204624237542,0.0155735319963907,0.152981516122118,0.114272108052331,0.391128516119853,0.33804157740185,4.11586803451685,5.95223619137436,0.0370925230615831,0.0877121507700289,0.0285364540447902,0.0256828086403112,-0.058608514682622,-0.0978214274139516,0.00702046242375134,0.0155735319963885,-1.39078799459603,-1.1142944367354
5,-0.342871767652598,-0.554506919015954,0.199940867464958,0.216657919742212,0.447147478428492,0.465465272326747,4.00014782860594,5.59764979346792,-0.0520194635641462,-0.0696089987726521,-0.202833831789998,-0.357877306170488,-0.0880184722984567,-0.127020614072812,-0.3428717676526,-0.554506919015952,-1.74068022467238,-1.68437488774775
6,0.191524645437639,0.440110207090828,0.188816172313118,0.217983462221048,0.434529829025716,0.466886990845802,4.02786252915087,5.59291172030795,0.00800973820069506,0.0500753419923305,0.115246263693708,0.274344181926019,0.0682686435432385,0.11569068317249,0.191524645437642,0.440110207090839,1.36441509216338,1.52523903712291
7,-0.308189084752417,-0.482966574707401,0.141437372423211,0.200445576913419,0.376081603409701,0.447711488476026,4.14382100862667,5.65527867778974,-0.099144837939897,-0.145827359673472,-0.208256444037528,-0.315253106749821,-0.000787802774987632,-0.0218861082841027,-0.308189084752413,-0.482966574707395,0.686828354441722,0
8,0.200678978955736,0.348856374480332,0.187187667563252,0.208313886948288,0.432651901143693,0.456414161643006,4.03190360112131,5.6273834525531,0.0249578986836898,0.0738802077728371,0.175721080272048,0.274976166707489,0,0,0.200678978955738,0.348856374480327,0,0
9,-0.018410346122461,-0.119367510725329,0.179508661013015,0.105647480996723,0.423684624470861,0.32503458430869,4.05090443971957,5.98114543880443,-0.112068758091788,-0.203660081497722,0.0936584119693256,0.0842925707723931,0,0,-0.0184103461224622,-0.119367510725329,0,0

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

저작자표시

'Animal Breeding > R for Genetic Evaluation' 카테고리의 다른 글

R을 이용한 일반적인 모형의 육종가 구하기 프로그래밍 (0)	2020.08.25
R을 이용한 다형질 개체 모형(unequal design matrices)의 육종가 구하기 프로그래밍 (0)	2020.07.29
R을 이용한 다형질 개체 모형(same model and no missing record)의 육종가 구하기 프로그래밍 (0)	2020.07.03
R을 이용하여 다형질 개체 모형(multiple trait animal model)에 대한 blupf90 분석 결과 정리하기 (0)	2020.07.02
R을 이용한 공통 환경 효과 모형(Common Environmental Effect Model)의 육종가 구하기 프로그래밍 (0)	2020.07.01

R을 이용한 다형질 개체 모형(same model and no missing record)의 육종가 구하기 프로그래밍

투정이 2020. 7. 3. 16:50

2020. 7. 3. 16:50

Multiple trait animal model, same model and no missing record

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

Raphael Mrode

Example 5.1 p72

자료

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

animal, sex, wwg(pre-weaning gain), pwg(post-weaning gain)

mt01_data.txt로 저장

혈통

자료와 혈통을 읽어 육종가를 구하는 R 프로그램이다.

신뢰도, 정확도, Standard error of prediction을 구한다.

육종가를 PA(parenta average), YD(yield devation), PC(progeny contribution)으로 구한다.

개체의 DYD(daughter yield deviation)을 구한다.

# multiple trait animal model, same model and no missing record - Date : 2020-07-03

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

# Raphael Mrode

# Example 5.1 p72

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

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

# functions

# make design matrix 
desgn <- function(v) {
    if (is.numeric(v)) {
        va = v
        mrow = length(va)
        mcol = max(va)
    }
    if (is.character(v)) {
        vf = factor(v)
        # 각 수준의 인덱스 값을 저장
        va = as.numeric(vf)
        mrow = length(va)
        mcol = length(levels(vf))
    }
    
    # 0으로 채워진 X 준비
    X = matrix(data = c(0), nrow = mrow, ncol = mcol)
    
    for (i in 1:mrow) {
        ic = va[i]
        X[i, ic] = 1
    }
    return(X)
}

# function to make inverse of numerator relationship matrix
ainv = function(pedi) {
    n = nrow(pedi)
    Ainv = matrix(c(0), nrow = n, ncol = n)
    
    for (i in 1:n) {
        animal = pedi[i, 1]
        sire = pedi[i, 2]
        dam = pedi[i, 3]
        
        if (sire == 0 & dam == 0) {
            # both parents unknown
            alpha = 1
            Ainv[animal, animal] = alpha + Ainv[animal, animal]
        } else if (sire != 0 & dam == 0) {
            # sire known
            alpha = 4/3
            Ainv[animal, animal] = alpha + Ainv[animal, animal]
            Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
            Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
            Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
        } else if (sire == 0 & dam != 0) {
            # dam known
            alpha = 4/3
            Ainv[animal, animal] = alpha + Ainv[animal, animal]
            Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
            Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
            Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
        } else {
            # both parents known
            alpha = 2
            Ainv[animal, animal] = alpha + Ainv[animal, animal]
            Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
            Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
            Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
            Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
            Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
            Ainv[sire, dam] = alpha/4 + Ainv[sire, dam]
            Ainv[dam, sire] = alpha/4 + Ainv[dam, sire]
            Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
        }
    }
    return(Ainv)
}

# set working directory 
setwd("D:/temp/04_multiple_traits_01_R") 

# print working directory 
getwd()

# no_of_trts
no_of_trts = 2

# prdigree and data

# read pedigree : animal sire dam
pedi = read.table("mt01_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("mt01_data.txt", header = FALSE, sep = "", col.names = c("animal", "sex", "wwg", "pwg"), stringsAsFactors = FALSE)

# print
data

# number of data
no_data = nrow(data)

# print
no_data

# variance component and ratio
G = matrix(c(20, 18, 18, 40), 2, 2)
R = matrix(c(40, 11, 11, 30), 2, 2)

# print
G
R

# inverse of G and R
Gi = ginv(G)
Ri = ginv(R)

# print
Gi
Ri

# full matrix of R inverse
Rif = kronecker(Ri, diag(no_data))

# print
Rif

# design matrix

# design matrix of 1st fixed effect
X1 = desgn(data[, 2])

# block diagonal matirx
X = as.matrix(bdiag(X1, X1))

# print
X

# number of levels of 1st fixed effect
no_lvl_f1 = ncol(X1)

# print
no_lvl_f1

# desing matrix of animal effect
Z1 = desgn(data[, 1])

# block diagonal matirx
Z = as.matrix(bdiag(Z1, Z1))

# print
Z

# number of animals
no_animals = ncol(Z1)

# print
no_animals

# observation
y1 = as.matrix(data[, 3])
y2 = as.matrix(data[, 4])
y = rbind(y1, y2)

# print
y1
y2
y

# inverse matrix of NRM
ai = ainv(pedi)

# print
ai

# LHS, RHS

# LHS construction
XPX = t(X) %*% Rif %*% X
XPZ = t(X) %*% Rif %*% Z
ZPX = t(Z) %*% Rif %*% X
ZPZ = t(Z) %*% Rif %*% Z

#print
XPX 
XPZ
ZPX
ZPZ

LHS = rbind(
        cbind(XPX, XPZ), 
        cbind(ZPX, ZPZ + kronecker(Gi, ai))
)

# print
LHS

# RHS construction
XPy = t(X) %*% Rif %*% y
ZPy = t(Z) %*% Rif %*% y

# print
XPy
ZPy

RHS = rbind(XPy,
            ZPy
)

# print
RHS

# Solutions

# generalized inverse of LHS
gi_LHS = ginv(LHS)

# print
gi_LHS

# solution
sol = gi_LHS %*% RHS

# print
sol

# solutions for fixed effects and animal effects
sol_t1_f1 = as.matrix(sol[1 : no_lvl_f1])
sol_t2_f1 = as.matrix(sol[(no_lvl_f1 + 1) : (no_lvl_f1 * 2)])
sol_f1 = cbind(sol_t1_f1, sol_t2_f1)

#print
sol_f1

# breedinv value
sol_t1_bv = sol[(no_lvl_f1 * 2 + 1) : (no_lvl_f1 * 2 + no_animals)]
sol_t2_bv = sol[(no_lvl_f1 * 2 + no_animals + 1) : (no_lvl_f1 * 2 + no_animals *2)]
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[(no_lvl_f1 * 2 + 1) : (no_lvl_f1 * 2 + no_animals), 
                       (no_lvl_f1 * 2 + 1) : (no_lvl_f1 * 2 + no_animals)])
d_ani_t2 = diag(gi_LHS[(no_lvl_f1 * 2 + no_animals + 1) : (no_lvl_f1 * 2 + no_animals * 2), 
                       (no_lvl_f1 * 2 + no_animals + 1) : (no_lvl_f1 * 2 + no_animals * 2)])
d_ani = cbind(d_ani_t1, d_ani_t2)

# print
d_ani

# reliability
rel = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)

for(i in 1 : no_animals){
  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, no_animals * no_of_trts), ncol = no_of_trts)

# numerator of n2
a2 = array(rep(0, no_of_trts * no_of_trts * no_animals), dim = c(no_of_trts, no_of_trts, no_animals))

# looping data
for (i in 1 : no_data) {
    yd1[data[i, 1], ] = as.matrix(data[i, c(3,4)] - sol_f1[data[i,2],])

    a2[,,data[i, 1]] = Ri
}
  
# print
yd1
a2

# Parents average, progeny contribution

# parents avearge
pa1 = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)

# progeny contribution numerator
pc0 = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)

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

# numerator of n1
a1 = array(rep(0, no_of_trts * no_of_trts * no_animals), dim = c(no_of_trts, no_of_trts, no_animals))


# looping pedi
for (i in 1 : no_animals) {

    sire = pedi[i, 2]
    dam = pedi[i, 3]
    
    if (sire == 0 & dam == 0) {
        # both parents unknown
        # PA
        a1[ , , i] = 1 * 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, no_animals * no_of_trts), ncol = no_of_trts)
for (i in 1 : no_animals) {
    pa[i, ] = ginv(n_de[ , , i]) %*% a1[, , i] %*% pa1[i, ]
}

# print
pa

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

# print
yd

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

# progeny contribution fraction of breeding values
pc = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
for (i in 1 : no_animals) {
  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_of_trts * no_of_trts * no_animals), dim = c(no_of_trts, no_of_trts, no_animals))
for (i in 1 : no_animals) {
  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, no_animals * no_of_trts), ncol = no_of_trts)
# denominator of dyd : summation of u of progeny * n2 of progeny
dyd_d = array(rep(0, no_of_trts * no_of_trts * no_animals), dim = c(no_of_trts, no_of_trts, no_animals))

# looping pedi
for (i in 1 : no_animals) {
# 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, no_animals * no_of_trts), ncol = no_of_trts)

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

# print
dyd

# breeding values and fractions
mt1_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
mt1_result

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

다음은 R 프로그램을 실행한 결과이다.

> # multiple trait animal model, same model and no missing record - Date : 2020-07-03
> 
> # Linear Models for the Prediction of Animal Breeding Values, 3rd Edtion
> 
> # Raphael Mrode
> 
> # Example 5.1 p72
> 
> # MASS packages 
> if (!("MASS" %in% installed.packages())){
+     install.packages('MASS', dependencies = TRUE)  
+ }
> library(MASS)
> 
> # Matrix packages 
> if (!("Matrix" %in% installed.packages())){
+   install.packages('Matrix', dependencies = TRUE)  
+ }
> library(Matrix)
> 
> # functions
> 
> # make design matrix 
> desgn <- function(v) {
+     if (is.numeric(v)) {
+         va = v
+         mrow = length(va)
+         mcol = max(va)
+     }
+     if (is.character(v)) {
+         vf = factor(v)
+         # 각 수준의 인덱스 값을 저장
+         va = as.numeric(vf)
+         mrow = length(va)
+         mcol = length(levels(vf))
+     }
+     
+     # 0으로 채워진 X 준비
+     X = matrix(data = c(0), nrow = mrow, ncol = mcol)
+     
+     for (i in 1:mrow) {
+         ic = va[i]
+         X[i, ic] = 1
+     }
+     return(X)
+ }
> 
> # function to make inverse of numerator relationship matrix
> ainv = function(pedi) {
+     n = nrow(pedi)
+     Ainv = matrix(c(0), nrow = n, ncol = n)
+     
+     for (i in 1:n) {
+         animal = pedi[i, 1]
+         sire = pedi[i, 2]
+         dam = pedi[i, 3]
+         
+         if (sire == 0 & dam == 0) {
+             # both parents unknown
+             alpha = 1
+             Ainv[animal, animal] = alpha + Ainv[animal, animal]
+         } else if (sire != 0 & dam == 0) {
+             # sire known
+             alpha = 4/3
+             Ainv[animal, animal] = alpha + Ainv[animal, animal]
+             Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
+             Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
+             Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
+         } else if (sire == 0 & dam != 0) {
+             # dam known
+             alpha = 4/3
+             Ainv[animal, animal] = alpha + Ainv[animal, animal]
+             Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
+             Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
+             Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
+         } else {
+             # both parents known
+             alpha = 2
+             Ainv[animal, animal] = alpha + Ainv[animal, animal]
+             Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
+             Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
+             Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
+             Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
+             Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
+             Ainv[sire, dam] = alpha/4 + Ainv[sire, dam]
+             Ainv[dam, sire] = alpha/4 + Ainv[dam, sire]
+             Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
+         }
+     }
+     return(Ainv)
+ }
> 
> # set working directory 
> setwd("D:/temp/04_multiple_traits_01_R") 
> 
> # print working directory 
> getwd()
[1] "D:/temp/04_multiple_traits_01_R"
> 
> # no_of_trts
> no_of_trts = 2
> 
> # prdigree and data
> 
> # read pedigree : animal sire dam
> pedi = read.table("mt01_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   0
5      5    3   2
6      6    1   2
7      7    4   5
8      8    3   6
> 
> # 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
> no_data = nrow(data)
> 
> # print
> no_data
[1] 5
> 
> # variance component and ratio
> G = matrix(c(20, 18, 18, 40), 2, 2)
> R = matrix(c(40, 11, 11, 30), 2, 2)
> 
> # print
> G
     [,1] [,2]
[1,]   20   18
[2,]   18   40
> R
     [,1] [,2]
[1,]   40   11
[2,]   11   30
> 
> # inverse of G and R
> Gi = ginv(G)
> Ri = ginv(R)
> 
> # print
> Gi
            [,1]        [,2]
[1,]  0.08403361 -0.03781513
[2,] -0.03781513  0.04201681
> Ri
            [,1]        [,2]
[1,]  0.02780352 -0.01019462
[2,] -0.01019462  0.03707136
> 
> # full matrix of R inverse
> Rif = kronecker(Ri, diag(no_data))
> 
> # print
> Rif
             [,1]        [,2]        [,3]        [,4]        [,5]        [,6]        [,7]
 [1,]  0.02780352  0.00000000  0.00000000  0.00000000  0.00000000 -0.01019462  0.00000000
 [2,]  0.00000000  0.02780352  0.00000000  0.00000000  0.00000000  0.00000000 -0.01019462
 [3,]  0.00000000  0.00000000  0.02780352  0.00000000  0.00000000  0.00000000  0.00000000
 [4,]  0.00000000  0.00000000  0.00000000  0.02780352  0.00000000  0.00000000  0.00000000
 [5,]  0.00000000  0.00000000  0.00000000  0.00000000  0.02780352  0.00000000  0.00000000
 [6,] -0.01019462  0.00000000  0.00000000  0.00000000  0.00000000  0.03707136  0.00000000
 [7,]  0.00000000 -0.01019462  0.00000000  0.00000000  0.00000000  0.00000000  0.03707136
 [8,]  0.00000000  0.00000000 -0.01019462  0.00000000  0.00000000  0.00000000  0.00000000
 [9,]  0.00000000  0.00000000  0.00000000 -0.01019462  0.00000000  0.00000000  0.00000000
[10,]  0.00000000  0.00000000  0.00000000  0.00000000 -0.01019462  0.00000000  0.00000000
             [,8]        [,9]       [,10]
 [1,]  0.00000000  0.00000000  0.00000000
 [2,]  0.00000000  0.00000000  0.00000000
 [3,] -0.01019462  0.00000000  0.00000000
 [4,]  0.00000000 -0.01019462  0.00000000
 [5,]  0.00000000  0.00000000 -0.01019462
 [6,]  0.00000000  0.00000000  0.00000000
 [7,]  0.00000000  0.00000000  0.00000000
 [8,]  0.03707136  0.00000000  0.00000000
 [9,]  0.00000000  0.03707136  0.00000000
[10,]  0.00000000  0.00000000  0.03707136
> 
> # design matrix
> 
> # design matrix of 1st fixed effect
> X1 = desgn(data[, 2])
> 
> # block diagonal matirx
> X = as.matrix(bdiag(X1, X1))
> 
> # print
> X
      [,1] [,2] [,3] [,4]
 [1,]    1    0    0    0
 [2,]    0    1    0    0
 [3,]    0    1    0    0
 [4,]    1    0    0    0
 [5,]    1    0    0    0
 [6,]    0    0    1    0
 [7,]    0    0    0    1
 [8,]    0    0    0    1
 [9,]    0    0    1    0
[10,]    0    0    1    0
> 
> # number of levels of 1st fixed effect
> no_lvl_f1 = ncol(X1)
> 
> # print
> no_lvl_f1
[1] 2
> 
> # desing matrix of animal effect
> Z1 = desgn(data[, 1])
> 
> # block diagonal matirx
> Z = as.matrix(bdiag(Z1, Z1))
> 
> # print
> Z
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
 [1,]    0    0    0    1    0    0    0    0    0     0     0     0     0     0     0     0
 [2,]    0    0    0    0    1    0    0    0    0     0     0     0     0     0     0     0
 [3,]    0    0    0    0    0    1    0    0    0     0     0     0     0     0     0     0
 [4,]    0    0    0    0    0    0    1    0    0     0     0     0     0     0     0     0
 [5,]    0    0    0    0    0    0    0    1    0     0     0     0     0     0     0     0
 [6,]    0    0    0    0    0    0    0    0    0     0     0     1     0     0     0     0
 [7,]    0    0    0    0    0    0    0    0    0     0     0     0     1     0     0     0
 [8,]    0    0    0    0    0    0    0    0    0     0     0     0     0     1     0     0
 [9,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     1     0
[10,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     1
> 
> # number of animals
> no_animals = ncol(Z1)
> 
> # print
> no_animals
[1] 8
> 
> # observation
> y1 = as.matrix(data[, 3])
> y2 = as.matrix(data[, 4])
> y = rbind(y1, y2)
> 
> # print
> y1
     [,1]
[1,]  4.5
[2,]  2.9
[3,]  3.9
[4,]  3.5
[5,]  5.0
> y2
     [,1]
[1,]  6.8
[2,]  5.0
[3,]  6.8
[4,]  6.0
[5,]  7.5
> y
      [,1]
 [1,]  4.5
 [2,]  2.9
 [3,]  3.9
 [4,]  3.5
 [5,]  5.0
 [6,]  6.8
 [7,]  5.0
 [8,]  6.8
 [9,]  6.0
[10,]  7.5
> 
> # inverse matrix of NRM
> ai = ainv(pedi)
> 
> # print
> ai
           [,1] [,2] [,3]       [,4] [,5] [,6] [,7] [,8]
[1,]  1.8333333  0.5  0.0 -0.6666667  0.0 -1.0    0    0
[2,]  0.5000000  2.0  0.5  0.0000000 -1.0 -1.0    0    0
[3,]  0.0000000  0.5  2.0  0.0000000 -1.0  0.5    0   -1
[4,] -0.6666667  0.0  0.0  1.8333333  0.5  0.0   -1    0
[5,]  0.0000000 -1.0 -1.0  0.5000000  2.5  0.0   -1    0
[6,] -1.0000000 -1.0  0.5  0.0000000  0.0  2.5    0   -1
[7,]  0.0000000  0.0  0.0 -1.0000000 -1.0  0.0    2    0
[8,]  0.0000000  0.0 -1.0  0.0000000  0.0 -1.0    0    2
> 
> # LHS, RHS
> 
> # LHS construction
> XPX = t(X) %*% Rif %*% X
> XPZ = t(X) %*% Rif %*% Z
> ZPX = t(Z) %*% Rif %*% X
> ZPZ = t(Z) %*% Rif %*% Z
> 
> #print
> XPX 
            [,1]        [,2]        [,3]        [,4]
[1,]  0.08341057  0.00000000 -0.03058387  0.00000000
[2,]  0.00000000  0.05560704  0.00000000 -0.02038925
[3,] -0.03058387  0.00000000  0.11121409  0.00000000
[4,]  0.00000000 -0.02038925  0.00000000  0.07414272
> XPZ
     [,1] [,2] [,3]        [,4]        [,5]        [,6]        [,7]        [,8] [,9] [,10]
[1,]    0    0    0  0.02780352  0.00000000  0.00000000  0.02780352  0.02780352    0     0
[2,]    0    0    0  0.00000000  0.02780352  0.02780352  0.00000000  0.00000000    0     0
[3,]    0    0    0 -0.01019462  0.00000000  0.00000000 -0.01019462 -0.01019462    0     0
[4,]    0    0    0  0.00000000 -0.01019462 -0.01019462  0.00000000  0.00000000    0     0
     [,11]       [,12]       [,13]       [,14]       [,15]       [,16]
[1,]     0 -0.01019462  0.00000000  0.00000000 -0.01019462 -0.01019462
[2,]     0  0.00000000 -0.01019462 -0.01019462  0.00000000  0.00000000
[3,]     0  0.03707136  0.00000000  0.00000000  0.03707136  0.03707136
[4,]     0  0.00000000  0.03707136  0.03707136  0.00000000  0.00000000
> ZPX
             [,1]        [,2]        [,3]        [,4]
 [1,]  0.00000000  0.00000000  0.00000000  0.00000000
 [2,]  0.00000000  0.00000000  0.00000000  0.00000000
 [3,]  0.00000000  0.00000000  0.00000000  0.00000000
 [4,]  0.02780352  0.00000000 -0.01019462  0.00000000
 [5,]  0.00000000  0.02780352  0.00000000 -0.01019462
 [6,]  0.00000000  0.02780352  0.00000000 -0.01019462
 [7,]  0.02780352  0.00000000 -0.01019462  0.00000000
 [8,]  0.02780352  0.00000000 -0.01019462  0.00000000
 [9,]  0.00000000  0.00000000  0.00000000  0.00000000
[10,]  0.00000000  0.00000000  0.00000000  0.00000000
[11,]  0.00000000  0.00000000  0.00000000  0.00000000
[12,] -0.01019462  0.00000000  0.03707136  0.00000000
[13,]  0.00000000 -0.01019462  0.00000000  0.03707136
[14,]  0.00000000 -0.01019462  0.00000000  0.03707136
[15,] -0.01019462  0.00000000  0.03707136  0.00000000
[16,] -0.01019462  0.00000000  0.03707136  0.00000000
> ZPZ
      [,1] [,2] [,3]        [,4]        [,5]        [,6]        [,7]        [,8] [,9] [,10]
 [1,]    0    0    0  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000    0     0
 [2,]    0    0    0  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000    0     0
 [3,]    0    0    0  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000    0     0
 [4,]    0    0    0  0.02780352  0.00000000  0.00000000  0.00000000  0.00000000    0     0
 [5,]    0    0    0  0.00000000  0.02780352  0.00000000  0.00000000  0.00000000    0     0
 [6,]    0    0    0  0.00000000  0.00000000  0.02780352  0.00000000  0.00000000    0     0
 [7,]    0    0    0  0.00000000  0.00000000  0.00000000  0.02780352  0.00000000    0     0
 [8,]    0    0    0  0.00000000  0.00000000  0.00000000  0.00000000  0.02780352    0     0
 [9,]    0    0    0  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000    0     0
[10,]    0    0    0  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000    0     0
[11,]    0    0    0  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000    0     0
[12,]    0    0    0 -0.01019462  0.00000000  0.00000000  0.00000000  0.00000000    0     0
[13,]    0    0    0  0.00000000 -0.01019462  0.00000000  0.00000000  0.00000000    0     0
[14,]    0    0    0  0.00000000  0.00000000 -0.01019462  0.00000000  0.00000000    0     0
[15,]    0    0    0  0.00000000  0.00000000  0.00000000 -0.01019462  0.00000000    0     0
[16,]    0    0    0  0.00000000  0.00000000  0.00000000  0.00000000 -0.01019462    0     0
      [,11]       [,12]       [,13]       [,14]       [,15]       [,16]
 [1,]     0  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000
 [2,]     0  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000
 [3,]     0  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000
 [4,]     0 -0.01019462  0.00000000  0.00000000  0.00000000  0.00000000
 [5,]     0  0.00000000 -0.01019462  0.00000000  0.00000000  0.00000000
 [6,]     0  0.00000000  0.00000000 -0.01019462  0.00000000  0.00000000
 [7,]     0  0.00000000  0.00000000  0.00000000 -0.01019462  0.00000000
 [8,]     0  0.00000000  0.00000000  0.00000000  0.00000000 -0.01019462
 [9,]     0  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000
[10,]     0  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000
[11,]     0  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000
[12,]     0  0.03707136  0.00000000  0.00000000  0.00000000  0.00000000
[13,]     0  0.00000000  0.03707136  0.00000000  0.00000000  0.00000000
[14,]     0  0.00000000  0.00000000  0.03707136  0.00000000  0.00000000
[15,]     0  0.00000000  0.00000000  0.00000000  0.03707136  0.00000000
[16,]     0  0.00000000  0.00000000  0.00000000  0.00000000  0.03707136
> 
> LHS = rbind(
+         cbind(XPX, XPZ), 
+         cbind(ZPX, ZPZ + kronecker(Gi, ai))
+ )
> 
> # print
> LHS
             [,1]        [,2]        [,3]        [,4]        [,5]        [,6]        [,7]
 [1,]  0.08341057  0.00000000 -0.03058387  0.00000000  0.00000000  0.00000000  0.00000000
 [2,]  0.00000000  0.05560704  0.00000000 -0.02038925  0.00000000  0.00000000  0.00000000
 [3,] -0.03058387  0.00000000  0.11121409  0.00000000  0.00000000  0.00000000  0.00000000
 [4,]  0.00000000 -0.02038925  0.00000000  0.07414272  0.00000000  0.00000000  0.00000000
 [5,]  0.00000000  0.00000000  0.00000000  0.00000000  0.15406162  0.04201681  0.00000000
 [6,]  0.00000000  0.00000000  0.00000000  0.00000000  0.04201681  0.16806723  0.04201681
 [7,]  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000  0.04201681  0.16806723
 [8,]  0.02780352  0.00000000 -0.01019462  0.00000000 -0.05602241  0.00000000  0.00000000
 [9,]  0.00000000  0.02780352  0.00000000 -0.01019462  0.00000000 -0.08403361 -0.08403361
[10,]  0.00000000  0.02780352  0.00000000 -0.01019462 -0.08403361 -0.08403361  0.04201681
[11,]  0.02780352  0.00000000 -0.01019462  0.00000000  0.00000000  0.00000000  0.00000000
[12,]  0.02780352  0.00000000 -0.01019462  0.00000000  0.00000000  0.00000000 -0.08403361
[13,]  0.00000000  0.00000000  0.00000000  0.00000000 -0.06932773 -0.01890756  0.00000000
[14,]  0.00000000  0.00000000  0.00000000  0.00000000 -0.01890756 -0.07563025 -0.01890756
[15,]  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000 -0.01890756 -0.07563025
[16,] -0.01019462  0.00000000  0.03707136  0.00000000  0.02521008  0.00000000  0.00000000
[17,]  0.00000000 -0.01019462  0.00000000  0.03707136  0.00000000  0.03781513  0.03781513
[18,]  0.00000000 -0.01019462  0.00000000  0.03707136  0.03781513  0.03781513 -0.01890756
[19,] -0.01019462  0.00000000  0.03707136  0.00000000  0.00000000  0.00000000  0.00000000
[20,] -0.01019462  0.00000000  0.03707136  0.00000000  0.00000000  0.00000000  0.03781513
             [,8]        [,9]       [,10]       [,11]       [,12]       [,13]       [,14]
 [1,]  0.02780352  0.00000000  0.00000000  0.02780352  0.02780352  0.00000000  0.00000000
 [2,]  0.00000000  0.02780352  0.02780352  0.00000000  0.00000000  0.00000000  0.00000000
 [3,] -0.01019462  0.00000000  0.00000000 -0.01019462 -0.01019462  0.00000000  0.00000000
 [4,]  0.00000000 -0.01019462 -0.01019462  0.00000000  0.00000000  0.00000000  0.00000000
 [5,] -0.05602241  0.00000000 -0.08403361  0.00000000  0.00000000 -0.06932773 -0.01890756
 [6,]  0.00000000 -0.08403361 -0.08403361  0.00000000  0.00000000 -0.01890756 -0.07563025
 [7,]  0.00000000 -0.08403361  0.04201681  0.00000000 -0.08403361  0.00000000 -0.01890756
 [8,]  0.18186515  0.04201681  0.00000000 -0.08403361  0.00000000  0.02521008  0.00000000
 [9,]  0.04201681  0.23788756  0.00000000 -0.08403361  0.00000000  0.00000000  0.03781513
[10,]  0.00000000  0.00000000  0.23788756  0.00000000 -0.08403361  0.03781513  0.03781513
[11,] -0.08403361 -0.08403361  0.00000000  0.19587075  0.00000000  0.00000000  0.00000000
[12,]  0.00000000  0.00000000 -0.08403361  0.00000000  0.19587075  0.00000000  0.00000000
[13,]  0.02521008  0.00000000  0.03781513  0.00000000  0.00000000  0.07703081  0.02100840
[14,]  0.00000000  0.03781513  0.03781513  0.00000000  0.00000000  0.02100840  0.08403361
[15,]  0.00000000  0.03781513 -0.01890756  0.00000000  0.03781513  0.00000000  0.02100840
[16,] -0.07952236 -0.01890756  0.00000000  0.03781513  0.00000000 -0.02801120  0.00000000
[17,] -0.01890756 -0.10473244  0.00000000  0.03781513  0.00000000  0.00000000 -0.04201681
[18,]  0.00000000  0.00000000 -0.10473244  0.00000000  0.03781513 -0.04201681 -0.04201681
[19,]  0.03781513  0.03781513  0.00000000 -0.08582488  0.00000000  0.00000000  0.00000000
[20,]  0.00000000  0.00000000  0.03781513  0.00000000 -0.08582488  0.00000000  0.00000000
            [,15]       [,16]       [,17]       [,18]       [,19]       [,20]
 [1,]  0.00000000 -0.01019462  0.00000000  0.00000000 -0.01019462 -0.01019462
 [2,]  0.00000000  0.00000000 -0.01019462 -0.01019462  0.00000000  0.00000000
 [3,]  0.00000000  0.03707136  0.00000000  0.00000000  0.03707136  0.03707136
 [4,]  0.00000000  0.00000000  0.03707136  0.03707136  0.00000000  0.00000000
 [5,]  0.00000000  0.02521008  0.00000000  0.03781513  0.00000000  0.00000000
 [6,] -0.01890756  0.00000000  0.03781513  0.03781513  0.00000000  0.00000000
 [7,] -0.07563025  0.00000000  0.03781513 -0.01890756  0.00000000  0.03781513
 [8,]  0.00000000 -0.07952236 -0.01890756  0.00000000  0.03781513  0.00000000
 [9,]  0.03781513 -0.01890756 -0.10473244  0.00000000  0.03781513  0.00000000
[10,] -0.01890756  0.00000000  0.00000000 -0.10473244  0.00000000  0.03781513
[11,]  0.00000000  0.03781513  0.03781513  0.00000000 -0.08582488  0.00000000
[12,]  0.03781513  0.00000000  0.00000000  0.03781513  0.00000000 -0.08582488
[13,]  0.00000000 -0.02801120  0.00000000 -0.04201681  0.00000000  0.00000000
[14,]  0.02100840  0.00000000 -0.04201681 -0.04201681  0.00000000  0.00000000
[15,]  0.08403361  0.00000000 -0.04201681  0.02100840  0.00000000 -0.04201681
[16,]  0.00000000  0.11410217  0.02100840  0.00000000 -0.04201681  0.00000000
[17,] -0.04201681  0.02100840  0.14211338  0.00000000 -0.04201681  0.00000000
[18,]  0.02100840  0.00000000  0.00000000  0.14211338  0.00000000 -0.04201681
[19,]  0.00000000 -0.04201681 -0.04201681  0.00000000  0.12110498  0.00000000
[20,] -0.04201681  0.00000000  0.00000000 -0.04201681  0.00000000  0.12110498
> 
> # RHS construction
> XPy = t(X) %*% Rif %*% y
> ZPy = t(Z) %*% Rif %*% y
> 
> # print
> XPy
           [,1]
[1,] 0.15449490
[2,] 0.06876738
[3,] 0.62001854
[4,] 0.36811863
> ZPy
            [,1]
 [1,] 0.00000000
 [2,] 0.00000000
 [3,] 0.00000000
 [4,] 0.05579240
 [5,] 0.02965709
 [6,] 0.03911029
 [7,] 0.03614458
 [8,] 0.06255792
 [9,] 0.00000000
[10,] 0.00000000
[11,] 0.00000000
[12,] 0.20620945
[13,] 0.15579240
[14,] 0.21232623
[15,] 0.18674699
[16,] 0.22706209
> 
> RHS = rbind(XPy,
+             ZPy
+ )
> 
> # print
> RHS
            [,1]
 [1,] 0.15449490
 [2,] 0.06876738
 [3,] 0.62001854
 [4,] 0.36811863
 [5,] 0.00000000
 [6,] 0.00000000
 [7,] 0.00000000
 [8,] 0.05579240
 [9,] 0.02965709
[10,] 0.03911029
[11,] 0.03614458
[12,] 0.06255792
[13,] 0.00000000
[14,] 0.00000000
[15,] 0.00000000
[16,] 0.20620945
[17,] 0.15579240
[18,] 0.21232623
[19,] 0.18674699
[20,] 0.22706209
> 
> # Solutions
> 
> # generalized inverse of LHS
> gi_LHS = ginv(LHS)
> 
> # print
> gi_LHS
            [,1]       [,2]       [,3]       [,4]       [,5]       [,6]       [,7]       [,8]
 [1,]  23.805131   6.291201  13.048276   5.671170 -6.5395332 -3.3374389 -5.2755138 -10.541038
 [2,]   6.291201  31.996882   5.671170  16.041950 -5.3702695 -9.5573818 -4.3877316  -3.712060
 [3,]  13.048276   5.671170  30.857577  12.600581 -5.8215059 -2.9950329 -4.8830435  -9.361379
 [4,]   5.671170  16.041950  12.600581  38.483277 -4.9970359 -8.3727339 -3.6798298  -3.924119
 [5,]  -6.539533  -5.370269  -5.821506  -4.997036 18.6047347  0.3243012  1.5791386   8.540106
 [6,]  -3.337439  -9.557382  -2.995033  -8.372734  0.3243012 19.5963751 -0.4802746   1.001548
 [7,]  -5.275514  -4.387732  -4.883044  -3.679830  1.5791386 -0.4802746 17.8928862   2.301944
 [8,] -10.541038  -3.712060  -9.361379  -3.924119  8.5401057  1.0015480  2.3019440  16.504656
 [9,]  -5.951896 -11.844286  -5.412983 -10.352380  2.1631118  9.3897226  7.6551424   2.269591
[10,]  -6.630505 -12.149478  -5.929358 -10.731520  8.5774271  9.7250410  1.1203208   5.154530
[11,] -11.291251  -7.376935  -9.968148  -6.529986  5.6324645  4.7003325  5.0739444   9.706120
[12,]  -9.583104  -7.784607  -8.815299  -6.559406  5.4460295  4.3104363  8.4506531   5.412337
[13,]  -5.821506  -4.997036 -12.950918 -11.068126 16.0904195  0.4332042  2.1649074   7.012413
[14,]  -2.995033  -8.372734  -6.657554 -18.656944  0.4332042 17.4241779 -0.6389804   1.388999
[15,]  -4.883044  -3.679830 -10.821190  -8.237206  2.1649074 -0.6389804 15.1103333   3.128063
[16,]  -9.361379  -3.924119 -20.830966  -8.590649  7.0124132  1.3889995  3.1280633  13.211993
[17,]  -5.412983 -10.352380 -12.016345 -23.073618  2.9533095  8.1651708  5.7955284   3.105642
[18,]  -5.929358 -10.731520 -13.184816 -23.892936  7.0407623  8.5802970  1.5641312   4.742595
[19,]  -9.968148  -6.529986 -22.194546 -14.535359  5.3386160  4.0369934  4.6333538   8.655863
[20,]  -8.815299  -6.559406 -19.547220 -14.675734  5.1134886  3.5591059  6.8877136   6.216282
            [,9]      [,10]      [,11]     [,12]       [,13]       [,14]       [,15]
 [1,]  -5.951896  -6.630505 -11.291251 -9.583104  -5.8215059  -2.9950329  -4.8830435
 [2,] -11.844286 -12.149478  -7.376935 -7.784607  -4.9970359  -8.3727339  -3.6798298
 [3,]  -5.412983  -5.929358  -9.968148 -8.815299 -12.9509185  -6.6575537 -10.8211898
 [4,] -10.352380 -10.731520  -6.529986 -6.559406 -11.0681257 -18.6569441  -8.2372059
 [5,]   2.163112   8.577427   5.632464  5.446029  16.0904195   0.4332042   2.1649074
 [6,]   9.389723   9.725041   4.700333  4.310436   0.4332042  17.4241779  -0.6389804
 [7,]   7.655142   1.120321   5.073944  8.450653   2.1649074  -0.6389804  15.1103333
 [8,]   2.269591   5.154530   9.706120  5.412337   7.0124132   1.3889995   3.1280633
 [9,]  16.541066   7.147506   8.548265  7.037833   2.9533095   8.1651708   5.7955284
[10,]   7.147506  17.151450   6.205604  8.531381   7.0407623   8.5802970   1.5641312
[11,]   8.548265   6.205604  17.115040  7.052592   5.3386160   4.0369934   4.6333538
[12,]   7.037833   8.531381   7.052592 16.284382   5.1134886   3.5591059   6.8877136
[13,]   2.953309   7.040762   5.338616  5.113489  35.9017859   0.9312686   4.6456425
[14,]   8.165171   8.580297   4.036993  3.559106   0.9312686  38.7676309  -1.3740157
[15,]   5.795528   1.564131   4.633354  6.887714   4.6456425  -1.3740157  33.7992438
[16,]   3.105642   4.742595   8.655863  6.216282  15.7328475   2.9783085   6.7165153
[17,]  13.278202   7.426558   7.024588  6.108719   6.3392413  18.2082861  13.1220853
[18,]   7.426558  14.036482   6.035384  7.010094  15.7970100  19.1056022   3.3523265
[19,]   7.024588   6.035384  13.970273  7.278308  11.8037249   9.0140533  10.2814964
[20,]   6.108719   7.010094   7.278308 12.951308  11.3161830   7.9802991  15.4655577
           [,16]      [,17]      [,18]      [,19]      [,20]
 [1,]  -9.361379  -5.412983  -5.929358  -9.968148  -8.815299
 [2,]  -3.924119 -10.352380 -10.731520  -6.529986  -6.559406
 [3,] -20.830966 -12.016345 -13.184816 -22.194546 -19.547220
 [4,]  -8.590649 -23.073618 -23.892936 -14.535359 -14.675734
 [5,]   7.012413   2.953309   7.040762   5.338616   5.113489
 [6,]   1.388999   8.165171   8.580297   4.036993   3.559106
 [7,]   3.128063   5.795528   1.564131   4.633354   6.887714
 [8,]  13.211993   3.105642   4.742595   8.655863   6.216282
 [9,]   3.105642  13.278202   7.426558   7.024588   6.108719
[10,]   4.742595   7.426558  14.036482   6.035384   7.010094
[11,]   8.655863   7.024588   6.035384  13.970273   7.278308
[12,]   6.216282   6.108719   7.010094   7.278308  12.951308
[13,]  15.732847   6.339241  15.797010  11.803725  11.316183
[14,]   2.978309  18.208286  19.105602   9.014053   7.980299
[15,]   6.716515  13.122085   3.352327  10.281496  15.465558
[16,]  29.724917   6.665202  10.516097  19.252951  13.515032
[17,]   6.665202  29.864618  16.282618  15.758829  13.625005
[18,]  10.516097  16.282618  31.503255  13.311889  15.726464
[19,]  19.252951  15.758829  13.311889  31.363554  15.967135
[20,]  13.515032  13.625005  15.726464  15.967135  29.159492
> 
> # solution
> sol = gi_LHS %*% RHS
> 
> # print
> sol
              [,1]
 [1,]  4.360866999
 [2,]  3.397261592
 [3,]  6.799897620
 [4,]  5.880295937
 [5,]  0.150915567
 [6,] -0.015392510
 [7,] -0.078391896
 [8,] -0.010238959
 [9,] -0.270331441
[10,]  0.275808258
[11,] -0.316117562
[12,]  0.243755523
[13,]  0.279597973
[14,] -0.007610071
[15,] -0.170341439
[16,] -0.012670709
[17,] -0.477830262
[18,]  0.517238387
[19,] -0.478983695
[20,]  0.391961543
> 
> # solutions for fixed effects and animal effects
> sol_t1_f1 = as.matrix(sol[1 : no_lvl_f1])
> sol_t2_f1 = as.matrix(sol[(no_lvl_f1 + 1) : (no_lvl_f1 * 2)])
> sol_f1 = cbind(sol_t1_f1, sol_t2_f1)
> 
> #print
> sol_f1
         [,1]     [,2]
[1,] 4.360867 6.799898
[2,] 3.397262 5.880296
> 
> # breedinv value
> sol_t1_bv = sol[(no_lvl_f1 * 2 + 1) : (no_lvl_f1 * 2 + no_animals)]
> sol_t2_bv = sol[(no_lvl_f1 * 2 + no_animals + 1) : (no_lvl_f1 * 2 + no_animals *2)]
> sol_bv = cbind(sol_t1_bv, sol_t2_bv)
> 
> # print
> sol_bv
       sol_t1_bv    sol_t2_bv
[1,]  0.15091557  0.279597973
[2,] -0.01539251 -0.007610071
[3,] -0.07839190 -0.170341439
[4,] -0.01023896 -0.012670709
[5,] -0.27033144 -0.477830262
[6,]  0.27580826  0.517238387
[7,] -0.31611756 -0.478983695
[8,]  0.24375552  0.391961543
> 
> # 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[(no_lvl_f1 * 2 + 1) : (no_lvl_f1 * 2 + no_animals), 
+                        (no_lvl_f1 * 2 + 1) : (no_lvl_f1 * 2 + no_animals)])
> d_ani_t2 = diag(gi_LHS[(no_lvl_f1 * 2 + no_animals + 1) : (no_lvl_f1 * 2 + no_animals * 2), 
+                        (no_lvl_f1 * 2 + no_animals + 1) : (no_lvl_f1 * 2 + no_animals * 2)])
> d_ani = cbind(d_ani_t1, d_ani_t2)
> 
> # print
> d_ani
     d_ani_t1 d_ani_t2
[1,] 18.60473 35.90179
[2,] 19.59638 38.76763
[3,] 17.89289 33.79924
[4,] 16.50466 29.72492
[5,] 16.54107 29.86462
[6,] 17.15145 31.50325
[7,] 17.11504 31.36355
[8,] 16.28438 29.15949
> 
> # reliability
> rel = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
> 
> for(i in 1 : no_animals){
+   rel[i, ] = 1 - d_ani[i,] / diag(G)  
+ }
> 
> # print
> rel
           [,1]       [,2]
[1,] 0.06976326 0.10245535
[2,] 0.02018124 0.03080923
[3,] 0.10535569 0.15501891
[4,] 0.17476719 0.25687708
[5,] 0.17294668 0.25338456
[6,] 0.14242750 0.21241863
[7,] 0.14424801 0.21591115
[8,] 0.18578089 0.27101269
> 
> # accuracy
> acc = sqrt(rel)
> 
> # print
> acc
          [,1]      [,2]
[1,] 0.2641274 0.3200865
[2,] 0.1420607 0.1755256
[3,] 0.3245854 0.3937244
[4,] 0.4180517 0.5068304
[5,] 0.4158686 0.5033732
[6,] 0.3773957 0.4608890
[7,] 0.3798000 0.4646624
[8,] 0.4310231 0.5205888
> 
> # standard error of prediction(SEP)
> sep = sqrt(d_ani)
> 
> # 확인
> sep
     d_ani_t1 d_ani_t2
[1,] 4.313321 5.991810
[2,] 4.426779 6.226366
[3,] 4.229998 5.813712
[4,] 4.062592 5.452056
[5,] 4.067071 5.464853
[6,] 4.141431 5.612776
[7,] 4.137033 5.600317
[8,] 4.035391 5.399953
> 
> # partitioning of breeding values
> 
> # yield deviation
> yd1 = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
> 
> # numerator of n2
> a2 = array(rep(0, no_of_trts * no_of_trts * no_animals), dim = c(no_of_trts, no_of_trts, no_animals))
> 
> # looping data
> for (i in 1 : no_data) {
+     yd1[data[i, 1], ] = as.matrix(data[i, c(3,4)] - sol_f1[data[i,2],])
+ 
+     a2[,,data[i, 1]] = Ri
+ }
>   
> # print
> yd1
           [,1]          [,2]
[1,]  0.0000000  0.0000000000
[2,]  0.0000000  0.0000000000
[3,]  0.0000000  0.0000000000
[4,]  0.1391330  0.0001023796
[5,] -0.4972616 -0.8802959373
[6,]  0.5027384  0.9197040627
[7,] -0.8608670 -0.7998976204
[8,]  0.6391330  0.7001023796
> 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.02780352 -0.01019462
[2,] -0.01019462  0.03707136

, , 5

            [,1]        [,2]
[1,]  0.02780352 -0.01019462
[2,] -0.01019462  0.03707136

, , 6

            [,1]        [,2]
[1,]  0.02780352 -0.01019462
[2,] -0.01019462  0.03707136

, , 7

            [,1]        [,2]
[1,]  0.02780352 -0.01019462
[2,] -0.01019462  0.03707136

, , 8

            [,1]        [,2]
[1,]  0.02780352 -0.01019462
[2,] -0.01019462  0.03707136

> 
> # Parents average, progeny contribution
> 
> # parents avearge
> pa1 = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
> 
> # progeny contribution numerator
> pc0 = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
> 
> # numerator of n3, denominator of progeny contribution
> a3 = array(rep(0, no_of_trts * no_of_trts * no_animals), dim = c(no_of_trts, no_of_trts, no_animals))
> 
> # numerator of n1
> a1 = array(rep(0, no_of_trts * no_of_trts * no_animals), dim = c(no_of_trts, no_of_trts, no_animals))
> 
> 
> # looping pedi
> for (i in 1 : no_animals) {
+ 
+     sire = pedi[i, 2]
+     dam = pedi[i, 3]
+     
+     if (sire == 0 & dam == 0) {
+         # both parents unknown
+         # PA
+         a1[ , , i] = 1 * 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.08403361 -0.03781513
[2,] -0.03781513  0.04201681

, , 2

            [,1]        [,2]
[1,]  0.08403361 -0.03781513
[2,] -0.03781513  0.04201681

, , 3

            [,1]        [,2]
[1,]  0.08403361 -0.03781513
[2,] -0.03781513  0.04201681

, , 4

            [,1]        [,2]
[1,]  0.11204482 -0.05042017
[2,] -0.05042017  0.05602241

, , 5

            [,1]        [,2]
[1,]  0.16806723 -0.07563025
[2,] -0.07563025  0.08403361

, , 6

            [,1]        [,2]
[1,]  0.16806723 -0.07563025
[2,] -0.07563025  0.08403361

, , 7

            [,1]        [,2]
[1,]  0.16806723 -0.07563025
[2,] -0.07563025  0.08403361

, , 8

            [,1]        [,2]
[1,]  0.16806723 -0.07563025
[2,] -0.07563025  0.08403361

> pa1
            [,1]        [,2]
[1,]  0.00000000  0.00000000
[2,]  0.00000000  0.00000000
[3,]  0.00000000  0.00000000
[4,]  0.07545778  0.13979899
[5,] -0.04689220 -0.08897575
[6,]  0.06776153  0.13599395
[7,] -0.14028520 -0.24525049
[8,]  0.09870818  0.17344847
> a3
, , 1

            [,1]        [,2]
[1,]  0.07002801 -0.03151261
[2,] -0.03151261  0.03501401

, , 2

            [,1]        [,2]
[1,]  0.08403361 -0.03781513
[2,] -0.03781513  0.04201681

, , 3

            [,1]        [,2]
[1,]  0.08403361 -0.03781513
[2,] -0.03781513  0.04201681

, , 4

            [,1]        [,2]
[1,]  0.04201681 -0.01890756
[2,] -0.01890756  0.02100840

, , 5

            [,1]        [,2]
[1,]  0.04201681 -0.01890756
[2,] -0.01890756  0.02100840

, , 6

            [,1]        [,2]
[1,]  0.04201681 -0.01890756
[2,] -0.01890756  0.02100840

, , 7

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

, , 8

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

> pc0
              [,1]          [,2]
[1,]  0.0038664044  0.0110750251
[2,] -0.0020114248  0.0005246376
[3,] -0.0002921427 -0.0083856077
[4,] -0.0061278140 -0.0032441978
[5,] -0.0082610361 -0.0080987423
[6,]  0.0057346179  0.0093477285
[7,]  0.0000000000  0.0000000000
[8,]  0.0000000000  0.0000000000
> 
> # denominator of n1, n2, n3, diagonal of animals in LHS
> n_de = a1 + a2 + a3
> 
> # print
> n_de
, , 1

            [,1]        [,2]
[1,]  0.15406162 -0.06932773
[2,] -0.06932773  0.07703081

, , 2

            [,1]        [,2]
[1,]  0.16806723 -0.07563025
[2,] -0.07563025  0.08403361

, , 3

            [,1]        [,2]
[1,]  0.16806723 -0.07563025
[2,] -0.07563025  0.08403361

, , 4

            [,1]        [,2]
[1,]  0.18186515 -0.07952236
[2,] -0.07952236  0.11410217

, , 5

           [,1]       [,2]
[1,]  0.2378876 -0.1047324
[2,] -0.1047324  0.1421134

, , 6

           [,1]       [,2]
[1,]  0.2378876 -0.1047324
[2,] -0.1047324  0.1421134

, , 7

            [,1]        [,2]
[1,]  0.19587075 -0.08582488
[2,] -0.08582488  0.12110498

, , 8

            [,1]        [,2]
[1,]  0.19587075 -0.08582488
[2,] -0.08582488  0.12110498

> 
> # parents average fraction of breeding values
> pa = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
> for (i in 1 : no_animals) {
+     pa[i, ] = ginv(n_de[ , , i]) %*% a1[, , i] %*% pa1[i, ]
+ }
> 
> # print
> pa
            [,1]        [,2]
[1,]  0.00000000  0.00000000
[2,]  0.00000000  0.00000000
[3,]  0.00000000  0.00000000
[4,]  0.03331733  0.05851558
[5,] -0.02519179 -0.04622283
[6,]  0.03577082  0.07071542
[7,] -0.08971192 -0.14614589
[8,]  0.06301831  0.10337078
> 
> # yield deviation fraction of breeding values
> yd = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
> for (i in 1 : no_animals) {
+   yd[i, ] = ginv(n_de[ , , i]) %*% a2[, , i] %*% yd1[i, ]
+ }
> 
> # print
> yd
           [,1]         [,2]
[1,]  0.0000000  0.000000000
[2,]  0.0000000  0.000000000
[3,]  0.0000000  0.000000000
[4,]  0.0227885  0.003484439
[5,] -0.1565945 -0.309364950
[6,]  0.1614856  0.322856538
[7,] -0.2264056 -0.332837809
[8,]  0.1807372  0.288590763
> 
> # progeny contribution
> pc1 = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
> for (i in 1 : no_animals) {
+   pc1[i, ] = ginv(a3[ , , i]) %*% pc0[i, ]
+ }
> pc1[is.nan(pc1) == TRUE] = 0
> pc1
            [,1]        [,2]
[1,]  0.33201425  0.61511554
[2,] -0.03078502 -0.01522014
[3,] -0.15678379 -0.34068288
[4,] -0.36190368 -0.48013713
[5,] -0.62199616 -0.94529668
[6,]  0.56590294  0.95426452
[7,]  0.00000000  0.00000000
[8,]  0.00000000  0.00000000
> 
> # progeny contribution fraction of breeding values
> pc = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
> for (i in 1 : no_animals) {
+   pc[i, ] = ginv(n_de[ , , i]) %*% a3[, , i] %*% pc1[i, ]
+ }
> 
> # print
> pc
            [,1]         [,2]
[1,]  0.15091557  0.279597973
[2,] -0.01539251 -0.007610071
[3,] -0.07839190 -0.170341439
[4,] -0.06634480 -0.074670727
[5,] -0.08854515 -0.122242479
[6,]  0.07855184  0.123666429
[7,]  0.00000000  0.000000000
[8,]  0.00000000  0.000000000
> 
> # Progeny(Daughter) Yield Deviation(PYD, DYD)
> 
> # n2 of progeny
> n2prog = array(rep(0, no_of_trts * no_of_trts * no_animals), dim = c(no_of_trts, no_of_trts, no_animals))
> for (i in 1 : no_animals) {
+   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.21085286 0.1389018
[2,] 0.02778035 0.4886564

, , 5

           [,1]      [,2]
[1,] 0.15238095 0.1190476
[2,] 0.02380952 0.3904762

, , 6

           [,1]      [,2]
[1,] 0.15238095 0.1190476
[2,] 0.02380952 0.3904762

, , 7

           [,1]      [,2]
[1,] 0.15238095 0.1190476
[2,] 0.02380952 0.3904762

, , 8

           [,1]      [,2]
[1,] 0.15238095 0.1190476
[2,] 0.02380952 0.3904762

> 
> # numerator of dyd : summation of u of progeny * n2 of progeny * (2 * YD - bv of mate)
> dyd_n = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
> # denominator of dyd : summation of u of progeny * n2 of progeny
> dyd_d = array(rep(0, no_of_trts * no_of_trts * no_animals), dim = c(no_of_trts, no_of_trts, no_animals))
> 
> # looping pedi
> for (i in 1 : no_animals) {
+ # 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,]  0.41457857  0.75074330
[2,] -0.01300594 -0.01335216
[3,] -0.10001866 -0.33916490
[4,] -0.35473336 -0.47265781
[5,] -0.44974264 -0.66048422
[6,]  0.39369861  0.64556227
[7,]  0.00000000  0.00000000
[8,]  0.00000000  0.00000000
> dyd_d
, , 1

           [,1]      [,2]
[1,] 0.29294952 0.2116488
[2,] 0.04232976 0.7162471

, , 2

           [,1]      [,2]
[1,] 0.30476190 0.2380952
[2,] 0.04761905 0.7809524

, , 3

           [,1]      [,2]
[1,] 0.30476190 0.2380952
[2,] 0.04761905 0.7809524

, , 4

           [,1]      [,2]
[1,] 0.15238095 0.1190476
[2,] 0.02380952 0.3904762

, , 5

           [,1]      [,2]
[1,] 0.15238095 0.1190476
[2,] 0.02380952 0.3904762

, , 6

           [,1]      [,2]
[1,] 0.15238095 0.1190476
[2,] 0.02380952 0.3904762

, , 7

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

, , 8

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

> 
> # dyd
> dyd = matrix(rep(0, no_animals * no_of_trts), ncol = no_of_trts)
> 
> # looping pedi
> for (i in 1 : no_animals) {
+     dyd[i, ] =  ginv(dyd_d[ , , i]) %*% dyd_n[i, ]
+ }  
> dyd[is.nan(dyd) == TRUE] = 0
> 
> # print
> dyd
            [,1]        [,2]
[1,]  0.68726086  1.00754574
[2,] -0.03078502 -0.01522014
[3,]  0.01166354 -0.43500772
[4,] -1.45140256 -1.12196498
[5,] -1.71149504 -1.58712453
[6,]  1.35665790  1.57054620
[7,]  0.00000000  0.00000000
[8,]  0.00000000  0.00000000
> 
> # breeding values and fractions
> mt1_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
> mt1_result
  animal animal_bv.sol_t1_bv animal_bv.sol_t2_bv      rel.1      rel.2     acc.1     acc.2
1      1          0.15091557         0.279597973 0.06976326 0.10245535 0.2641274 0.3200865
2      2         -0.01539251        -0.007610071 0.02018124 0.03080923 0.1420607 0.1755256
3      3         -0.07839190        -0.170341439 0.10535569 0.15501891 0.3245854 0.3937244
4      4         -0.01023896        -0.012670709 0.17476719 0.25687708 0.4180517 0.5068304
5      5         -0.27033144        -0.477830262 0.17294668 0.25338456 0.4158686 0.5033732
6      6          0.27580826         0.517238387 0.14242750 0.21241863 0.3773957 0.4608890
7      7         -0.31611756        -0.478983695 0.14424801 0.21591115 0.3798000 0.4646624
8      8          0.24375552         0.391961543 0.18578089 0.27101269 0.4310231 0.5205888
  sep.d_ani_t1 sep.d_ani_t2        pa.1        pa.2       yd.1         yd.2        pc.1
1     4.313321     5.991810  0.00000000  0.00000000  0.0000000  0.000000000  0.15091557
2     4.426779     6.226366  0.00000000  0.00000000  0.0000000  0.000000000 -0.01539251
3     4.229998     5.813712  0.00000000  0.00000000  0.0000000  0.000000000 -0.07839190
4     4.062592     5.452056  0.03331733  0.05851558  0.0227885  0.003484439 -0.06634480
5     4.067071     5.464853 -0.02519179 -0.04622283 -0.1565945 -0.309364950 -0.08854515
6     4.141431     5.612776  0.03577082  0.07071542  0.1614856  0.322856538  0.07855184
7     4.137033     5.600317 -0.08971192 -0.14614589 -0.2264056 -0.332837809  0.00000000
8     4.035391     5.399953  0.06301831  0.10337078  0.1807372  0.288590763  0.00000000
          pc.2 sum_of_fr.1  sum_of_fr.2       dyd.1       dyd.2
1  0.279597973  0.15091557  0.279597973  0.68726086  1.00754574
2 -0.007610071 -0.01539251 -0.007610071 -0.03078502 -0.01522014
3 -0.170341439 -0.07839190 -0.170341439  0.01166354 -0.43500772
4 -0.074670727 -0.01023896 -0.012670709 -1.45140256 -1.12196498
5 -0.122242479 -0.27033144 -0.477830262 -1.71149504 -1.58712453
6  0.123666429  0.27580826  0.517238387  1.35665790  1.57054620
7  0.000000000 -0.31611756 -0.478983695  0.00000000  0.00000000
8  0.000000000  0.24375552  0.391961543  0.00000000  0.00000000
> 
> # 파일로 출력, 분리자는 ",", 따옴표 없고 
> output_filename = gsub("[ -]", "", paste("mt1_result_", Sys.Date(), ".csv")) 
> write.table(mt1_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,0.150915567312797,0.279597972744865,0.0697632637735797,0.102455353618367,0.264127362788447,0.320086478343536,4.31332061462261,5.99180989812472,0,0,0,0,0.150915567312795,0.279597972744862,0.150915567312795,0.279597972744862,0.68726085711682,1.00754574375344
2,-0.0153925096914798,-0.0076100708199065,0.0201812427467539,0.030809227182284,0.14206070092307,0.175525574154549,4.42677931967078,6.22636578693451,0,0,0,0,-0.0153925096914804,-0.00761007081990604,-0.0153925096914804,-0.00761007081990604,-0.0307850193829532,-0.015220141639821
3,-0.0783918961641927,-0.170341438593948,0.105355689244771,0.155018905719144,0.324585411324617,0.39372440325581,4.22999837057942,5.81371170348464,0,0,0,0,-0.078391896164192,-0.170341438593946,-0.078391896164192,-0.170341438593946,0.0116635350966923,-0.435007715776464
4,-0.0102389585292884,-0.0126707086241291,0.174767187326832,0.256877083095507,0.418051656290024,0.506830428344143,4.06259230706988,5.45205618791477,0.0333173344095044,0.0585155786066336,0.0227885035278452,0.00348443941480283,-0.0663447964666395,-0.0746707266455604,-0.0102389585292899,-0.0126707086241239,-1.45140255671963,-1.12196497916736
5,-0.270331441389235,-0.477830261602733,0.172946680277051,0.253384555790556,0.415868585345239,0.50337317746435,4.06707098468406,5.46485295029773,-0.0251917915640488,-0.0462228319310023,-0.156594500241952,-0.3093649501734,-0.0885451495832349,-0.122242479498329,-0.270331441389235,-0.477830261602731,-1.71149503957958,-1.58712453214596
6,0.275808257580577,0.517238387038379,0.142427497991221,0.212418627668353,0.377395678289009,0.460888953727851,4.14143091698698,5.61277604160952,0.0357708240803181,0.070715419872911,0.161485595241754,0.322856537923106,0.078551838258502,0.123666429242362,0.275808257580574,0.517238387038379,1.35665789805533,1.57054619782385
7,-0.316117561639437,-0.478983695055088,0.144248005041,0.215911154973304,0.379799953977091,0.46466240968396,4.13703274088809,5.60031729467785,-0.0897119212614896,-0.146145886165347,-0.226405640377943,-0.332837808889743,0,0,-0.316117561639433,-0.47898369505509,0,0
8,0.2437555230054,0.391961542524077,0.185780891703132,0.271012693707549,0.431023075604001,0.52058879521898,4.03539120358081,5.39995298606368,0.0630183062404894,0.103370779985251,0.180737216764914,0.288590762538829,0,0,0.243755523005404,0.391961542524079,0,0

엑셀에서 열면 다음과 같다.

저작자표시

'Animal Breeding > R for Genetic Evaluation' 카테고리의 다른 글

R을 이용한 다형질 개체 모형(unequal design matrices)의 육종가 구하기 프로그래밍 (0)	2020.07.29
R을 이용한 다형질 개체 모형(same model and missing record)의 육종가 구하기 프로그래밍 (0)	2020.07.26
R을 이용하여 다형질 개체 모형(multiple trait animal model)에 대한 blupf90 분석 결과 정리하기 (0)	2020.07.02
R을 이용한 공통 환경 효과 모형(Common Environmental Effect Model)의 육종가 구하기 프로그래밍 (0)	2020.07.01
R을 이용하여 공통 환경 효과 모형(common environmental effect model)에 대한 blupf90 분석 결과 정리하기 (0)	2020.07.01

R을 이용하여 다형질 개체 모형(multiple trait animal model)에 대한 blupf90 분석 결과 정리하기

투정이 2020. 7. 2. 17:24

2020. 7. 2. 17:24

blupf90을 이용하여 다형질 개체 모형(다변량 개체 모형)의 육종가를 구하는 것을 아래에서 설명한 바 있다.

2014/02/25 - [Animal Breeding/BLUPF90] - blupf90으로 다형질 모형(same model and no missing record)의 육종가 구하기

2014/03/12 - [Animal Breeding/BLUPF90] - blupf90으로 다형질 모형(same model with missing record)의 육종가 구하기

2014/03/12 - [Animal Breeding/BLUPF90] - blupf90으로 다형질 모형(different model)의 육종가 구하기

2014/03/13 - [Animal Breeding/BLUPF90] - blupf90으로 다형질 모형(no environmental covariance)의 육종가 구하기

여기서는 blupf90이 구한 결과를 R로 정리하는 것을 설명한다.

blupf90은 육종가(breeding value)와 육종가의 SE(SEP : Standard Error of Prediction)를 제공하는데 SE를 이용하여 정확도(accuracy)를 구해야 한다. 또한 renumber된 개체 ID이기 때문에 원래 ID를 찾아 주어야 한다.

blupf90이 구한 solutions는 다음과 같다.

trait/effect level  solution          s.e.
   1   1         1          4.36086700          4.87905017
   2   1         1          6.79989768          5.55495970
   1   1         2          3.39726169          5.65657865
   2   1         2          5.88029605          6.20348910
   1   2         1         -0.31611757          4.13703273
   2   2         1         -0.47898372          5.60031728
   1   2         2         -0.01023894          4.06259227
   2   2         2         -0.01267064          5.45205614
   1   2         3          0.27580827          4.14143092
   2   2         3          0.51723842          5.61277604
   1   2         4          0.24375552          4.03539120
   2   2         4          0.39196151          5.39995298
   1   2         5         -0.27033146          4.06707098
   2   2         5         -0.47783034          5.46485295
   1   2         6          0.15091558          4.31332061
   2   2         6          0.27959802          5.99180990
   1   2         7         -0.01539252          4.42677932
   2   2         7         -0.00761009          6.22636579
   1   2         8         -0.07839191          4.22999837
   2   2         8         -0.17034149          5.81371170

renumber된 혈통과 원래 ID는 renadd02.ped에 나와 있다.

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

1, 2, 3열이 리넘버된 혈통이고 10열이 원래 ID이다.

다형질 개체 모형(multiple trait animal model)에 대한 blupf90의 결과로부터 정확도를 계산하고 원래 ID를 붙이는 R 프로그램이다.

# multiple traits animal model의 blupf90의 분석 결과를 R로 정리하기 - Date :2020-07-02 

# blupf90의 육종가와 SE(standard error of prediction)읽어 

#정확도를 계산하고, 원래 ID를 붙임 

# set working directory 

setwd("D:/temp/04_multiple_traits_01_blupf90") 

# print working directory 

getwd() 

# sigma_a 

sigma_a_t1 = 20 
sigma_a_t2 = 40 

sigma_a_t1
sigma_a_t2

# read solutions 

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

solutions 

# read renumbered pedigree 

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

pedi 

# read renumbered data 

data = read.table("renf90.dat", sep = "", stringsAsFactors = FALSE) 

data 

# read fixed solutions from solutions 

sol_fixed_t1 = solutions[(solutions$trait == 1 & solutions$effect == 1), 4] 

sol_fixed_t1

sol_fixed_t2 = solutions[(solutions$trait == 2 & solutions$effect == 1), 4] 

sol_fixed_t2

# read animal solutions from solutions 

sol_animal_t1 = solutions[(solutions$trait == 1 & solutions$effect == 2), 3:5] 

rownames(sol_animal_t1) = sol_animal_t1[,1] 

sol_animal_t1

sol_animal_t2 = solutions[(solutions$trait == 2 & solutions$effect == 2), 3:5] 

rownames(sol_animal_t2) = sol_animal_t2[,1] 

sol_animal_t2

# reliability 

rel_t1 = 1 - (sol_animal_t1$se^2 / sigma_a_t1) 

rel_t1

rel_t2 = 1 - (sol_animal_t2$se^2 / sigma_a_t2) 

rel_t2

# accuracy 

acc_t1 = sqrt(rel_t1) 

acc_t1 

acc_t2 = sqrt(rel_t2) 

acc_t2

# 육종가 분할하기, DYD 구하기 생략 

# result 

mt01_result = data.frame(animal = sol_animal_t1$level, bv_t1 = sol_animal_t1$solution, rel_t1 = rel_t1, acc_t1 = acc_t1, sep_t1 = sol_animal_t1$se,
                                                       bv_t2 = sol_animal_t2$solution, rel_t2 = rel_t2, acc_t2 = acc_t2, sep_t2 = sol_animal_t2$se)


mt01_result 

# join original animal id 

aid_renum = pedi[,c(1,10)] 

colnames(aid_renum) = c("animal", "orig_aid") 

aid_renum 

# 조인하기 위한 패키지 실행 

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

mt01_final = join(x = mt01_result, y = aid_renum, by = "animal", type = "left") 

mt01_final = mt01_final[order(mt01_final$orig_aid),] 

mt01_final

# 파일로 출력, 분리자는 ",", 따옴표 없고 

output_filename = gsub("[ -]", "", paste("mt01_result_", Sys.Date(), ".csv")) 

write.table(mt01_final, output_filename, sep=",", row.names = FALSE, quote = FALSE)

다음은 실행 결과이다.

> # multiple traits animal model의 blupf90의 분석 결과를 R로 정리하기 - Date :2020-07-02 
> 
> # blupf90의 육종가와 SE(standard error of prediction)읽어 
> 
> #정확도를 계산하고, 원래 ID를 붙임 
> 
> # set working directory 
> 
> setwd("D:/temp/04_multiple_traits_01_blupf90") 
> 
> # print working directory 
> 
> getwd() 
[1] "D:/temp/04_multiple_traits_01_blupf90"
> 
> # sigma_a 
> 
> sigma_a_t1 = 20 
> sigma_a_t2 = 40 
> 
> sigma_a_t1
[1] 20
> sigma_a_t2
[1] 40
> 
> # read solutions 
> 
> solutions = read.table("solutions", skip = 1, sep = "", col.names = c("trait", "effect", "level", "solution", "se"), stringsAsFactors = FALSE) 
> 
> solutions 
   trait effect level    solution       se
1      1      1     1  4.36086700 4.879050
2      2      1     1  6.79989768 5.554960
3      1      1     2  3.39726169 5.656579
4      2      1     2  5.88029605 6.203489
5      1      2     1 -0.31611757 4.137033
6      2      2     1 -0.47898372 5.600317
7      1      2     2 -0.01023894 4.062592
8      2      2     2 -0.01267064 5.452056
9      1      2     3  0.27580827 4.141431
10     2      2     3  0.51723842 5.612776
11     1      2     4  0.24375552 4.035391
12     2      2     4  0.39196151 5.399953
13     1      2     5 -0.27033146 4.067071
14     2      2     5 -0.47783034 5.464853
15     1      2     6  0.15091558 4.313321
16     2      2     6  0.27959802 5.991810
17     1      2     7 -0.01539252 4.426779
18     2      2     7 -0.00761009 6.226366
19     1      2     8 -0.07839191 4.229998
20     2      2     8 -0.17034149 5.813712
> 
> # 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
> 
> # read renumbered data 
> 
> data = read.table("renf90.dat", sep = "", stringsAsFactors = FALSE) 
> 
> data 
   V1  V2 V3 V4
1 4.5 6.8  1  2
2 2.9 5.0  2  5
3 3.9 6.8  2  3
4 3.5 6.0  1  1
5 5.0 7.5  1  4
> 
> # read fixed solutions from solutions 
> 
> sol_fixed_t1 = solutions[(solutions$trait == 1 & solutions$effect == 1), 4] 
> 
> sol_fixed_t1
[1] 4.360867 3.397262
> 
> sol_fixed_t2 = solutions[(solutions$trait == 2 & solutions$effect == 1), 4] 
> 
> sol_fixed_t2
[1] 6.799898 5.880296
> 
> # read animal solutions from solutions 
> 
> sol_animal_t1 = solutions[(solutions$trait == 1 & solutions$effect == 2), 3:5] 
> 
> rownames(sol_animal_t1) = sol_animal_t1[,1] 
> 
> sol_animal_t1
  level    solution       se
1     1 -0.31611757 4.137033
2     2 -0.01023894 4.062592
3     3  0.27580827 4.141431
4     4  0.24375552 4.035391
5     5 -0.27033146 4.067071
6     6  0.15091558 4.313321
7     7 -0.01539252 4.426779
8     8 -0.07839191 4.229998
> 
> sol_animal_t2 = solutions[(solutions$trait == 2 & solutions$effect == 2), 3:5] 
> 
> rownames(sol_animal_t2) = sol_animal_t2[,1] 
> 
> sol_animal_t2
  level    solution       se
1     1 -0.47898372 5.600317
2     2 -0.01267064 5.452056
3     3  0.51723842 5.612776
4     4  0.39196151 5.399953
5     5 -0.47783034 5.464853
6     6  0.27959802 5.991810
7     7 -0.00761009 6.226366
8     8 -0.17034149 5.813712
> 
> # reliability 
> 
> rel_t1 = 1 - (sol_animal_t1$se^2 / sigma_a_t1) 
> 
> rel_t1
[1] 0.14424801 0.17476720 0.14242750 0.18578089 0.17294668 0.06976327 0.02018124 0.10535569
> 
> rel_t2 = 1 - (sol_animal_t2$se^2 / sigma_a_t2) 
> 
> rel_t2
[1] 0.21591116 0.25687710 0.21241863 0.27101270 0.25338456 0.10245535 0.03080923 0.15501891
> 
> # accuracy 
> 
> acc_t1 = sqrt(rel_t1) 
> 
> acc_t1 
[1] 0.3798000 0.4180517 0.3773957 0.4310231 0.4158686 0.2641274 0.1420607 0.3245854
> 
> acc_t2 = sqrt(rel_t2) 
> 
> acc_t2
[1] 0.4646624 0.5068304 0.4608890 0.5205888 0.5033732 0.3200865 0.1755256 0.3937244
> 
> # 육종가 분할하기, DYD 구하기 생략 
> 
> # result 
> 
> mt01_result = data.frame(animal = sol_animal_t1$level, bv_t1 = sol_animal_t1$solution, rel_t1 = rel_t1, acc_t1 = acc_t1, sep_t1 = sol_animal_t1$se,
+                                                        bv_t2 = sol_animal_t2$solution, rel_t2 = rel_t2, acc_t2 = acc_t2, sep_t2 = sol_animal_t2$se)
> 
> 
> mt01_result 
  animal       bv_t1     rel_t1    acc_t1   sep_t1       bv_t2     rel_t2    acc_t2   sep_t2
1      1 -0.31611757 0.14424801 0.3798000 4.137033 -0.47898372 0.21591116 0.4646624 5.600317
2      2 -0.01023894 0.17476720 0.4180517 4.062592 -0.01267064 0.25687710 0.5068304 5.452056
3      3  0.27580827 0.14242750 0.3773957 4.141431  0.51723842 0.21241863 0.4608890 5.612776
4      4  0.24375552 0.18578089 0.4310231 4.035391  0.39196151 0.27101270 0.5205888 5.399953
5      5 -0.27033146 0.17294668 0.4158686 4.067071 -0.47783034 0.25338456 0.5033732 5.464853
6      6  0.15091558 0.06976327 0.2641274 4.313321  0.27959802 0.10245535 0.3200865 5.991810
7      7 -0.01539252 0.02018124 0.1420607 4.426779 -0.00761009 0.03080923 0.1755256 6.226366
8      8 -0.07839191 0.10535569 0.3245854 4.229998 -0.17034149 0.15501891 0.3937244 5.813712
> 
> # join original animal id 
> 
> 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
> 
> # 조인하기 위한 패키지 실행 
> 
> #install.packages("plyr") 
> if (!("plyr" %in% installed.packages())){
+   install.packages('plyr', dependencies = TRUE)  
+ }
> library(plyr) 
> 
> mt01_final = join(x = mt01_result, y = aid_renum, by = "animal", type = "left") 
> 
> mt01_final = mt01_final[order(mt01_final$orig_aid),] 
> 
> mt01_final
  animal       bv_t1     rel_t1    acc_t1   sep_t1       bv_t2     rel_t2    acc_t2   sep_t2 orig_aid
6      6  0.15091558 0.06976327 0.2641274 4.313321  0.27959802 0.10245535 0.3200865 5.991810        1
7      7 -0.01539252 0.02018124 0.1420607 4.426779 -0.00761009 0.03080923 0.1755256 6.226366        2
8      8 -0.07839191 0.10535569 0.3245854 4.229998 -0.17034149 0.15501891 0.3937244 5.813712        3
2      2 -0.01023894 0.17476720 0.4180517 4.062592 -0.01267064 0.25687710 0.5068304 5.452056        4
5      5 -0.27033146 0.17294668 0.4158686 4.067071 -0.47783034 0.25338456 0.5033732 5.464853        5
3      3  0.27580827 0.14242750 0.3773957 4.141431  0.51723842 0.21241863 0.4608890 5.612776        6
1      1 -0.31611757 0.14424801 0.3798000 4.137033 -0.47898372 0.21591116 0.4646624 5.600317        7
4      4  0.24375552 0.18578089 0.4310231 4.035391  0.39196151 0.27101270 0.5205888 5.399953        8
> 
> # 파일로 출력, 분리자는 ",", 따옴표 없고 
> 
> output_filename = gsub("[ -]", "", paste("mt01_result_", Sys.Date(), ".csv")) 
> 
> write.table(mt01_final, output_filename, sep=",", row.names = FALSE, quote = FALSE)

다음은 저장한 파일이다.

animal,bv_t1,rel_t1,acc_t1,sep_t1,bv_t2,rel_t2,acc_t2,sep_t2,orig_aid
6,0.15091558,0.0697632657674614,0.264127366562917,4.31332061,0.27959802,0.10245535305655,0.320086477465934,5.9918099,1
7,-0.01539252,0.0201812426010169,0.142060700410131,4.42677932,-0.00761009,0.030809226227942,0.175525571436022,6.22636579,2
8,-0.07839191,0.105355689489867,0.32458541170217,4.22999837,-0.17034149,0.155018906732078,0.393724404542159,5.8137117,3
2,-0.01023894,0.174767202386813,0.418051674302128,4.06259227,-0.01267064,0.256877096157208,0.506830441229814,5.45205614,4
5,-0.27033146,0.172946682182092,0.415868587635676,4.06707098,-0.47783034,0.253384555871907,0.503373177545156,5.46485295,5
3,0.27580827,0.142427496743397,0.377395676635806,4.14143092,0.51723842,0.212418628120048,0.460888954217877,5.61277604,6
1,-0.31611757,0.144248009545437,0.379799959907104,4.13703273,-0.47898372,0.215911159083335,0.464662414106559,5.60031728,7
4,0.24375552,0.185780893148128,0.43102307728024,4.0353912,0.39196151,0.271012695344728,0.52058879679141,5.39995298,8

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

저작자표시

'Animal Breeding > R for Genetic Evaluation' 카테고리의 다른 글

R을 이용한 다형질 개체 모형(same model and missing record)의 육종가 구하기 프로그래밍 (0)	2020.07.26
R을 이용한 다형질 개체 모형(same model and no missing record)의 육종가 구하기 프로그래밍 (0)	2020.07.03
R을 이용한 공통 환경 효과 모형(Common Environmental Effect Model)의 육종가 구하기 프로그래밍 (0)	2020.07.01
R을 이용하여 공통 환경 효과 모형(common environmental effect model)에 대한 blupf90 분석 결과 정리하기 (0)	2020.07.01
R을 이용한 반복 모형(Repeatability Model)의 육종가 구하기 프로그래밍 (0)	2020.06.14

PREV 이전 1 ···4 5 6 7 8 9 10 ···50 NEXT 다음

Let's Study Animal Breeding

전체 글

pgAdmin 4를 이용한 Backup과 Restore

'ETC' 카테고리의 다른 글

PostgreSQL psql을 이용한 role, database, table 생성

'ETC' 카테고리의 다른 글

R을 이용한 다형질 개체 모형(same model and missing record)의 육종가 구하기 프로그래밍

'Animal Breeding > R for Genetic Evaluation' 카테고리의 다른 글

R을 이용한 다형질 개체 모형(same model and no missing record)의 육종가 구하기 프로그래밍

'Animal Breeding > R for Genetic Evaluation' 카테고리의 다른 글

R을 이용하여 다형질 개체 모형(multiple trait animal model)에 대한 blupf90 분석 결과 정리하기

'Animal Breeding > R for Genetic Evaluation' 카테고리의 다른 글

+ Recent posts

티스토리툴바