# Single-step GBLUP Model(ssGBLUP)
# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.
# Raphael Mrode
# Example 11.6 p192
single step GBLUP(ssGBLUP)은 유전체 자료를 가지고 있는 개체와 가지고 있지 않은 개체를 동시에 분석할 수 있는 방법이다. 또한 SNP efffect를 구한 다음 DGV 또는 GEBV를 구하는 것이 아니라 SNP effect를 구하는 단계를 생략하고 한 번에 유전체 자료를 가지고 있는 개체는 GEBV, 없는 개체는 EBV를 구한다.
ssGBLUP에 포함하는 개체는 다음의 세 가지로 분류할 수 있다.
1. 능력검정 기록만 가지고 있는 개체 - 유전체 자료를 생성하기 전에, 즉 과거에 유전평가에 포함된 개체
2. 능력검정 기록과 유전체 자료를 가지고 있는 개체 - 소위 참조집단(reference population)이라고 하는 개체들로 GEBV를 추정한다. 참조집단이 많아야 유전체 선발(genomic selection)의 효율(아래 3번 개체들의 정확도 향상)이 높아진다.
3. 유전체 자료만 가지고 있는 개체 - 보통 어린 개체들로 아직 능력검정을 할 단계에 이르지 못한 개체. 이들에 대해 GEBV를 계산하고 조기 선발함으로써 세대간격(generation interval)을 줄여 유전적 개량량을 높일 수 있다. 이들은 시간이 지나 능력검정할 단계에 다다르면 능력검정을 하여 표현형 자료를 확보하여 참조집단이 된다.
모형 또는 MME는 위의 책 또는 아래 포스팅을 참고한다.
masuday.github.io/blupf90_tutorial/mrode_c11ex116_ssgblup.html
Numerical examples from Mrode (2014)
Numerical examples from Mrode (2014) Yutaka Masuda September 2019 Back to index.html. Single-step approach Model We consider the standard animal model \[ \mathbf{y}=\mathbf{Xb}+\mathbf{Wu}+\mathbf{e}. \] If some animals are genotyped, their additive relati
masuday.github.io
Data(ssgblup_data.txt)
animal sire dam mean edc fat_dyd snp1 snp2 snp3 snp4 snp5 snp6 snp7 snp8 snp9 snp10
13 0 0 1 558 9 2 0 1 1 0 0 0 2 1 2
14 0 0 1 722 13.4 1 0 0 0 0 2 0 2 1 0
15 13 4 1 300 12.7 1 1 2 1 1 0 0 2 1 2
16 15 2 1 73 15.4 0 0 2 1 0 1 0 2 2 1
17 15 5 1 52 5.9 0 1 1 2 0 0 0 2 1 2
18 14 6 1 87 7.7 1 1 0 1 0 2 0 2 2 1
19 14 9 1 64 10.2 0 0 1 1 0 2 0 2 2 0
20 14 9 1 103 4.8 0 1 1 0 0 1 0 2 2 0
21 1 3 1 13 7.6 2 0 0 0 0 1 2 2 1 2
22 14 8 1 125 8.8 0 0 0 1 1 2 0 2 0 0
23 14 11 1 93 9.8 0 1 1 0 0 1 0 2 2 1
24 14 10 1 66 9.2 1 0 0 0 1 1 0 2 0 0
25 14 7 1 75 11.5 0 0 0 1 1 2 0 2 1 0
26 14 12 1 33 13.3 1 0 1 1 0 2 0 1 0 0
Pedigree(ssgblup_pedi.txt)
animal sire dam
1 0 0
2 0 0
3 0 0
4 0 0
5 0 0
6 0 0
7 0 0
8 0 0
9 0 0
10 0 0
11 0 0
12 0 0
13 0 0
14 0 0
15 13 4
16 15 2
17 15 5
18 14 6
19 14 9
20 14 9
21 1 3
22 14 8
23 14 11
24 14 10
25 14 7
26 14 12
R 프로그램
# single step GBLUP Model
# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.
# Raphael Mrode
# Example 11.6 p192
# MASS packages
if (!("MASS" %in% installed.packages())) {
install.packages("MASS", dependencies = TRUE)
}
library(MASS)
# 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)
}
# function to make numerator relationship matrix
makenrm = function(pedi) {
n = nrow(pedi)
nrm = 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
nrm[animal, animal] = 1
} else if (sire != 0 & dam == 0) {
# sire known
for (j in 1:animal - 1) {
nrm[j, animal] = 0.5 * (nrm[j, sire])
nrm[animal, j] = nrm[j, animal]
}
nrm[animal, animal] = 1
} else if (sire == 0 & dam != 0) {
# dam known
for (j in 1:animal - 1) {
nrm[j, animal] = 0.5 * (nrm[j, dam])
nrm[animal, j] = nrm[j, animal]
}
nrm[animal, animal] = 1
} else {
# both parents known
for (j in 1:animal - 1) {
nrm[j, animal] = 0.5 * (nrm[j, sire] + nrm[j, dam])
nrm[animal, j] = nrm[j, animal]
}
nrm[animal, animal] = 1 + 0.5 * nrm[sire, dam]
}
}
return(nrm)
}
# set working directory
setwd(".")
# print working directory
getwd()
# read data
data = read.table("ssgblup_data.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
data
# read pedigree : animal sire dam
pedi = read.table("ssgblup_pedi.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
pedi
# variance component and ratio
sigma_a = 35.241
sigma_e = 245
sigma_a
sigma_e
# variance ratio, alpha
alpha = sigma_e / sigma_a
alpha
# SNP
M_all = as.matrix(data[6:14, 7:16])
M_all
# mean of each SNP(allele frequency)
snp_mean = colMeans(M_all) / 2
snp_mean
# sum of 2pq(heterozygote ratio)
sum2pq = sum(2 * snp_mean * (1 - snp_mean))
sum2pq
# observation
y = data[1:10, 6]
y
# design matrix of fixed effect(general mean)
X = matrix(rep(1, 10), 10, 1)
X
# design matrix for animal(parents : 12, animals with record : 10, animals without record : 4)
W = cbind(matrix(c(0), 10, 12), diag(10), matrix(c(0), 10, 4))
W
# genomic relationship matrix, G
Z = sweep(M_all, 2, snp_mean * 2)
round(Z, 3)
G = Z %*% t(Z) / sum2pq
round(G, 3)
# Numerator Relationship Matrix
A = makenrm(pedi)
A
A22 = A[18:26, 18:26]
A22
# weighted G
wt = 0.95
Gw = wt * G + (1 - wt) * A22
round(Gw, 3)
Gwi = ginv(Gw)
round(Gwi, 3)
# inverse matrix of NRM
Ai = ainv(pedi)
Ai
# H matrix
H0 = matrix(c(0), 26, 26)
H0[18:26, 18:26] = Gwi - ginv(A22)
H0
Hi = Ai + H0
round(Hi, 3)
# weight R
R = diag(data[1:10, 5])
R
Ri = ginv(R)
Ri
# Left Hand Side
LHS = rbind(
cbind(t(X) %*% Ri %*% X, t(X) %*% Ri %*% W),
cbind(t(W) %*% Ri %*% X, t(W) %*% Ri %*% W + Hi * alpha))
round(LHS, 3)
# Right Hand Side
RHS = rbind(t(X) %*% Ri %*% y,
t(W) %*% Ri %*% y)
round(RHS, 3)
# Solutions
# generalized inverse of LHS
gi_LHS = ginv(LHS)
round(gi_LHS, 3)
# solution
sol = gi_LHS %*% RHS
round(sol, 6)
# general mean effect
gme = sol[1]
gme
# (g)ebv of animals
ebv = sol[2:27]
ebv
RStudio에서 실행하는 화면
실행 결과
> # single step GBLUP Model
>
> # Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.
>
> # Raphael Mrode
>
> # Example 11.6 p192
>
> # MASS packages
> if (!("MASS" %in% installed.packages())) {
+ install.packages("MASS", dependencies = TRUE)
+ }
> library(MASS)
>
> # 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)
+ }
>
> # function to make numerator relationship matrix
> makenrm = function(pedi) {
+ n = nrow(pedi)
+ nrm = 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
+ nrm[animal, animal] = 1
+ } else if (sire != 0 & dam == 0) {
+ # sire known
+ for (j in 1:animal - 1) {
+ nrm[j, animal] = 0.5 * (nrm[j, sire])
+ nrm[animal, j] = nrm[j, animal]
+ }
+ nrm[animal, animal] = 1
+ } else if (sire == 0 & dam != 0) {
+ # dam known
+ for (j in 1:animal - 1) {
+ nrm[j, animal] = 0.5 * (nrm[j, dam])
+ nrm[animal, j] = nrm[j, animal]
+ }
+ nrm[animal, animal] = 1
+ } else {
+ # both parents known
+ for (j in 1:animal - 1) {
+ nrm[j, animal] = 0.5 * (nrm[j, sire] + nrm[j, dam])
+ nrm[animal, j] = nrm[j, animal]
+ }
+ nrm[animal, animal] = 1 + 0.5 * nrm[sire, dam]
+ }
+ }
+ return(nrm)
+ }
>
> # set working directory
> setwd(".")
>
> # print working directory
> getwd()
[1] "D:/users/bhpark/2021/job/공부하기/01_Linear Models for the Prediction of Animal Breeding Values_3rd_Mrode/25_ssGBLUP"
>
> # read data
> data = read.table("ssgblup_data.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
> data
animal sire dam mean edc fat_dyd snp1 snp2 snp3 snp4 snp5 snp6 snp7 snp8 snp9 snp10
1 13 0 0 1 558 9.0 2 0 1 1 0 0 0 2 1 2
2 14 0 0 1 722 13.4 1 0 0 0 0 2 0 2 1 0
3 15 13 4 1 300 12.7 1 1 2 1 1 0 0 2 1 2
4 16 15 2 1 73 15.4 0 0 2 1 0 1 0 2 2 1
5 17 15 5 1 52 5.9 0 1 1 2 0 0 0 2 1 2
6 18 14 6 1 87 7.7 1 1 0 1 0 2 0 2 2 1
7 19 14 9 1 64 10.2 0 0 1 1 0 2 0 2 2 0
8 20 14 9 1 103 4.8 0 1 1 0 0 1 0 2 2 0
9 21 1 3 1 13 7.6 2 0 0 0 0 1 2 2 1 2
10 22 14 8 1 125 8.8 0 0 0 1 1 2 0 2 0 0
11 23 14 11 1 93 9.8 0 1 1 0 0 1 0 2 2 1
12 24 14 10 1 66 9.2 1 0 0 0 1 1 0 2 0 0
13 25 14 7 1 75 11.5 0 0 0 1 1 2 0 2 1 0
14 26 14 12 1 33 13.3 1 0 1 1 0 2 0 1 0 0
>
> # read pedigree : animal sire dam
> pedi = read.table("ssgblup_pedi.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
> pedi
animal sire dam
1 1 0 0
2 2 0 0
3 3 0 0
4 4 0 0
5 5 0 0
6 6 0 0
7 7 0 0
8 8 0 0
9 9 0 0
10 10 0 0
11 11 0 0
12 12 0 0
13 13 0 0
14 14 0 0
15 15 13 4
16 16 15 2
17 17 15 5
18 18 14 6
19 19 14 9
20 20 14 9
21 21 1 3
22 22 14 8
23 23 14 11
24 24 14 10
25 25 14 7
26 26 14 12
>
> # variance component and ratio
> sigma_a = 35.241
> sigma_e = 245
> sigma_a
[1] 35.241
> sigma_e
[1] 245
>
> # variance ratio, alpha
> alpha = sigma_e / sigma_a
> alpha
[1] 6.95213
>
> # SNP
> M_all = as.matrix(data[6:14, 7:16])
> M_all
snp1 snp2 snp3 snp4 snp5 snp6 snp7 snp8 snp9 snp10
6 1 1 0 1 0 2 0 2 2 1
7 0 0 1 1 0 2 0 2 2 0
8 0 1 1 0 0 1 0 2 2 0
9 2 0 0 0 0 1 2 2 1 2
10 0 0 0 1 1 2 0 2 0 0
11 0 1 1 0 0 1 0 2 2 1
12 1 0 0 0 1 1 0 2 0 0
13 0 0 0 1 1 2 0 2 1 0
14 1 0 1 1 0 2 0 1 0 0
>
> # mean of each SNP(allele frequency)
> snp_mean = colMeans(M_all) / 2
> snp_mean
snp1 snp2 snp3 snp4 snp5 snp6 snp7 snp8 snp9 snp10
0.2777778 0.1666667 0.2222222 0.2777778 0.1666667 0.7777778 0.1111111 0.9444444 0.5555556 0.2222222
>
> # sum of 2pq(heterozygote ratio)
> sum2pq = sum(2 * snp_mean * (1 - snp_mean))
> sum2pq
[1] 3.191358
>
> # observation
> y = data[1:10, 6]
> y
[1] 9.0 13.4 12.7 15.4 5.9 7.7 10.2 4.8 7.6 8.8
>
> # design matrix of fixed effect(general mean)
> X = matrix(rep(1, 10), 10, 1)
> X
[,1]
[1,] 1
[2,] 1
[3,] 1
[4,] 1
[5,] 1
[6,] 1
[7,] 1
[8,] 1
[9,] 1
[10,] 1
>
> # design matrix for animal(parents : 12, animals with record : 10, animals without record : 4)
> W = cbind(matrix(c(0), 10, 12), diag(10), matrix(c(0), 10, 4))
> W
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
[1,] 0 0 0 0 0 0 0 0 0 0 0 0 1 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 1 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 1 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 1 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 1 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 1 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 1 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 1 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 1 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 1 0 0 0 0
>
> # genomic relationship matrix, G
> Z = sweep(M_all, 2, snp_mean * 2)
> round(Z, 3)
snp1 snp2 snp3 snp4 snp5 snp6 snp7 snp8 snp9 snp10
6 0.444 0.667 -0.444 0.444 -0.333 0.444 -0.222 0.111 0.889 0.556
7 -0.556 -0.333 0.556 0.444 -0.333 0.444 -0.222 0.111 0.889 -0.444
8 -0.556 0.667 0.556 -0.556 -0.333 -0.556 -0.222 0.111 0.889 -0.444
9 1.444 -0.333 -0.444 -0.556 -0.333 -0.556 1.778 0.111 -0.111 1.556
10 -0.556 -0.333 -0.444 0.444 0.667 0.444 -0.222 0.111 -1.111 -0.444
11 -0.556 0.667 0.556 -0.556 -0.333 -0.556 -0.222 0.111 0.889 0.556
12 0.444 -0.333 -0.444 -0.556 0.667 -0.556 -0.222 0.111 -1.111 -0.444
13 -0.556 -0.333 -0.444 0.444 0.667 0.444 -0.222 0.111 -0.111 -0.444
14 0.444 -0.333 0.556 0.444 -0.333 0.444 -0.222 -0.889 -1.111 -0.444
> G = Z %*% t(Z) / sum2pq
> round(G, 3)
6 7 8 9 10 11 12 13 14
6 0.785 0.124 0.054 0.193 -0.398 0.228 -0.538 -0.120 -0.329
7 0.124 0.716 0.333 -0.781 -0.120 0.193 -0.573 0.159 -0.050
8 0.054 0.333 0.890 -0.538 -0.503 0.750 -0.329 -0.224 -0.433
9 0.193 -0.781 -0.538 2.735 -0.677 -0.050 0.124 -0.712 -0.294
10 -0.398 -0.120 -0.503 -0.677 0.925 -0.642 0.472 0.576 0.368
11 0.228 0.193 0.750 -0.050 -0.642 0.925 -0.468 -0.364 -0.573
12 -0.538 -0.573 -0.329 0.124 0.472 -0.468 0.959 0.124 0.228
13 -0.120 0.159 -0.224 -0.712 0.576 -0.364 0.124 0.542 0.019
14 -0.329 -0.050 -0.433 -0.294 0.368 -0.573 0.228 0.019 1.064
>
> # Numerator Relationship Matrix
> A = makenrm(pedi)
> A
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
[1,] 1.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00 0.0 0.0 0.00 0.00 0.00 0.00 0.00 0.5 0.00 0.00 0.00 0.00 0.00
[2,] 0.0 1.0 0.0 0.00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00 0.0 0.0 0.50 0.00 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00
[3,] 0.0 0.0 1.0 0.00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00 0.0 0.0 0.00 0.00 0.00 0.00 0.00 0.5 0.00 0.00 0.00 0.00 0.00
[4,] 0.0 0.0 0.0 1.00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00 0.0 0.5 0.25 0.25 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00
[5,] 0.0 0.0 0.0 0.00 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00 0.0 0.0 0.00 0.50 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00
[6,] 0.0 0.0 0.0 0.00 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00 0.0 0.0 0.00 0.00 0.50 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00
[7,] 0.0 0.0 0.0 0.00 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.00 0.0 0.0 0.00 0.00 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.50 0.00
[8,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.00 0.0 0.0 0.00 0.00 0.00 0.00 0.00 0.0 0.50 0.00 0.00 0.00 0.00
[9,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.00 0.0 0.0 0.00 0.00 0.00 0.50 0.50 0.0 0.00 0.00 0.00 0.00 0.00
[10,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.00 0.0 0.0 0.00 0.00 0.00 0.00 0.00 0.0 0.00 0.00 0.50 0.00 0.00
[11,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.00 0.0 0.0 0.00 0.00 0.00 0.00 0.00 0.0 0.00 0.50 0.00 0.00 0.00
[12,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.00 0.0 0.0 0.00 0.00 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.50
[13,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.00 0.0 0.5 0.25 0.25 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00
[14,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00 1.0 0.0 0.00 0.00 0.50 0.50 0.50 0.0 0.50 0.50 0.50 0.50 0.50
[15,] 0.0 0.0 0.0 0.50 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.50 0.0 1.0 0.50 0.50 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00
[16,] 0.0 0.5 0.0 0.25 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.25 0.0 0.5 1.00 0.25 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00
[17,] 0.0 0.0 0.0 0.25 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.25 0.0 0.5 0.25 1.00 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00
[18,] 0.0 0.0 0.0 0.00 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.0 0.00 0.5 0.0 0.00 0.00 1.00 0.25 0.25 0.0 0.25 0.25 0.25 0.25 0.25
[19,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.00 0.5 0.0 0.00 0.00 0.25 1.00 0.50 0.0 0.25 0.25 0.25 0.25 0.25
[20,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.00 0.5 0.0 0.00 0.00 0.25 0.50 1.00 0.0 0.25 0.25 0.25 0.25 0.25
[21,] 0.5 0.0 0.5 0.00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00 0.0 0.0 0.00 0.00 0.00 0.00 0.00 1.0 0.00 0.00 0.00 0.00 0.00
[22,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0 0.00 0.5 0.0 0.00 0.00 0.25 0.25 0.25 0.0 1.00 0.25 0.25 0.25 0.25
[23,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.00 0.5 0.0 0.00 0.00 0.25 0.25 0.25 0.0 0.25 1.00 0.25 0.25 0.25
[24,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.00 0.5 0.0 0.00 0.00 0.25 0.25 0.25 0.0 0.25 0.25 1.00 0.25 0.25
[25,] 0.0 0.0 0.0 0.00 0.0 0.0 0.5 0.0 0.0 0.0 0.0 0.0 0.00 0.5 0.0 0.00 0.00 0.25 0.25 0.25 0.0 0.25 0.25 0.25 1.00 0.25
[26,] 0.0 0.0 0.0 0.00 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.5 0.00 0.5 0.0 0.00 0.00 0.25 0.25 0.25 0.0 0.25 0.25 0.25 0.25 1.00
> A22 = A[18:26, 18:26]
> A22
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 1.00 0.25 0.25 0 0.25 0.25 0.25 0.25 0.25
[2,] 0.25 1.00 0.50 0 0.25 0.25 0.25 0.25 0.25
[3,] 0.25 0.50 1.00 0 0.25 0.25 0.25 0.25 0.25
[4,] 0.00 0.00 0.00 1 0.00 0.00 0.00 0.00 0.00
[5,] 0.25 0.25 0.25 0 1.00 0.25 0.25 0.25 0.25
[6,] 0.25 0.25 0.25 0 0.25 1.00 0.25 0.25 0.25
[7,] 0.25 0.25 0.25 0 0.25 0.25 1.00 0.25 0.25
[8,] 0.25 0.25 0.25 0 0.25 0.25 0.25 1.00 0.25
[9,] 0.25 0.25 0.25 0 0.25 0.25 0.25 0.25 1.00
>
> # weighted G
> wt = 0.95
> Gw = wt * G + (1 - wt) * A22
> round(Gw, 3)
6 7 8 9 10 11 12 13 14
6 0.796 0.130 0.064 0.184 -0.366 0.229 -0.498 -0.101 -0.300
7 0.130 0.730 0.341 -0.742 -0.101 0.196 -0.531 0.163 -0.035
8 0.064 0.341 0.895 -0.511 -0.465 0.725 -0.300 -0.201 -0.399
9 0.184 -0.742 -0.511 2.648 -0.643 -0.048 0.118 -0.676 -0.279
10 -0.366 -0.101 -0.465 -0.643 0.928 -0.598 0.461 0.560 0.362
11 0.229 0.196 0.725 -0.048 -0.598 0.928 -0.432 -0.333 -0.531
12 -0.498 -0.531 -0.300 0.118 0.461 -0.432 0.961 0.130 0.229
13 -0.101 0.163 -0.201 -0.676 0.560 -0.333 0.130 0.565 0.031
14 -0.300 -0.035 -0.399 -0.279 0.362 -0.531 0.229 0.031 1.061
> Gwi = ginv(Gw)
> round(Gwi, 3)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 2.535 0.934 0.907 0.410 1.522 0.065 1.368 -0.812 0.438
[2,] 0.934 5.634 -1.782 0.567 1.995 1.049 3.148 -3.458 -0.805
[3,] 0.907 -1.782 6.846 1.825 2.193 -3.183 -1.582 1.542 1.208
[4,] 0.410 0.567 1.825 1.593 1.227 0.217 0.018 1.324 0.889
[5,] 1.522 1.995 2.193 1.227 7.496 -0.183 -0.323 -5.482 -0.773
[6,] 0.065 1.049 -3.183 0.217 -0.183 5.225 1.512 1.685 1.217
[7,] 1.368 3.148 -1.582 0.018 -0.323 1.512 3.965 -0.903 -0.062
[8,] -0.812 -3.458 1.542 1.324 -5.482 1.685 -0.903 11.227 3.166
[9,] 0.438 -0.805 1.208 0.889 -0.773 1.217 -0.062 3.166 2.523
>
> # inverse matrix of NRM
> Ai = ainv(pedi)
> Ai
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
[1,] 1.5 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.0 0 0 0 0 0 -1 0 0 0 0 0
[2,] 0.0 1.5 0.0 0.0 0.0 0.0 0.0 0.0 0 0.0 0.0 0.0 0.0 0.0 0.5 -1 0 0 0 0 0 0 0 0 0 0
[3,] 0.5 0.0 1.5 0.0 0.0 0.0 0.0 0.0 0 0.0 0.0 0.0 0.0 0.0 0.0 0 0 0 0 0 -1 0 0 0 0 0
[4,] 0.0 0.0 0.0 1.5 0.0 0.0 0.0 0.0 0 0.0 0.0 0.0 0.5 0.0 -1.0 0 0 0 0 0 0 0 0 0 0 0
[5,] 0.0 0.0 0.0 0.0 1.5 0.0 0.0 0.0 0 0.0 0.0 0.0 0.0 0.0 0.5 0 -1 0 0 0 0 0 0 0 0 0
[6,] 0.0 0.0 0.0 0.0 0.0 1.5 0.0 0.0 0 0.0 0.0 0.0 0.0 0.5 0.0 0 0 -1 0 0 0 0 0 0 0 0
[7,] 0.0 0.0 0.0 0.0 0.0 0.0 1.5 0.0 0 0.0 0.0 0.0 0.0 0.5 0.0 0 0 0 0 0 0 0 0 0 -1 0
[8,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.5 0 0.0 0.0 0.0 0.0 0.5 0.0 0 0 0 0 0 0 -1 0 0 0 0
[9,] 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 1.0 0.0 0 0 0 -1 -1 0 0 0 0 0 0
[10,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0 1.5 0.0 0.0 0.0 0.5 0.0 0 0 0 0 0 0 0 0 -1 0 0
[11,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0 0.0 1.5 0.0 0.0 0.5 0.0 0 0 0 0 0 0 0 -1 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 1.5 0.0 0.5 0.0 0 0 0 0 0 0 0 0 0 0 -1
[13,] 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0 0 0.0 0.0 0.0 1.5 0.0 -1.0 0 0 0 0 0 0 0 0 0 0 0
[14,] 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.5 1 0.5 0.5 0.5 0.0 5.0 0.0 0 0 -1 -1 -1 0 -1 -1 -1 -1 -1
[15,] 0.0 0.5 0.0 -1.0 0.5 0.0 0.0 0.0 0 0.0 0.0 0.0 -1.0 0.0 3.0 -1 -1 0 0 0 0 0 0 0 0 0
[16,] 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0 0.0 0.0 0.0 0.0 0.0 -1.0 2 0 0 0 0 0 0 0 0 0 0
[17,] 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0 0.0 0.0 0.0 0.0 0.0 -1.0 0 2 0 0 0 0 0 0 0 0 0
[18,] 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0 0.0 0.0 0.0 0.0 -1.0 0.0 0 0 2 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 -1 0.0 0.0 0.0 0.0 -1.0 0.0 0 0 0 2 0 0 0 0 0 0 0
[20,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1 0.0 0.0 0.0 0.0 -1.0 0.0 0 0 0 0 2 0 0 0 0 0 0
[21,] -1.0 0.0 -1.0 0.0 0.0 0.0 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
[22,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0 0.0 0.0 0.0 0.0 -1.0 0.0 0 0 0 0 0 0 2 0 0 0 0
[23,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0 0.0 -1.0 0.0 0.0 -1.0 0.0 0 0 0 0 0 0 0 2 0 0 0
[24,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0 -1.0 0.0 0.0 0.0 -1.0 0.0 0 0 0 0 0 0 0 0 2 0 0
[25,] 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0 0.0 0.0 0.0 0.0 -1.0 0.0 0 0 0 0 0 0 0 0 0 2 0
[26,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0 0.0 0.0 -1.0 0.0 -1.0 0.0 0 0 0 0 0 0 0 0 0 0 2
>
> # H matrix
> H0 = matrix(c(0), 26, 26)
> H0[18:26, 18:26] = Gwi - ginv(A22)
> H0
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
[1,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000 0.000000 0.00000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[2,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000 0.000000 0.00000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[3,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000 0.000000 0.00000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[4,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000 0.000000 0.00000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[5,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000 0.000000 0.00000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[6,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000 0.000000 0.00000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[7,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000 0.000000 0.00000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[8,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000 0.000000 0.00000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[9,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000 0.000000 0.00000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[10,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000 0.000000 0.00000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[11,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000 0.000000 0.00000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[12,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000 0.000000 0.00000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[13,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000 0.000000 0.00000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[14,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000 0.000000 0.00000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[15,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000 0.000000 0.00000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[16,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000 0.000000 0.00000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[17,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0000000 0.0000000 0.000000 0.00000000 0.0000000 0.0000000 0.00000000 0.0000000 0.00000000
[18,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.3290432 1.0290474 1.002121 0.40953260 1.6487596 0.1916196 1.49505551 -0.6854842 0.56535657
[19,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0290474 4.2055366 -1.210784 0.56677591 2.0907120 1.1447303 3.24326156 -3.3624887 -0.70927508
[20,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0021206 -1.2107838 5.417277 1.82466493 2.2884902 -3.0879351 -1.48680074 1.6374007 1.30339943
[21,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.4095326 0.5667759 1.824665 0.59303100 1.2270063 0.2169856 0.01800482 1.3244139 0.88865550
[22,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.6487596 2.0907120 2.288490 1.22700625 6.2894776 -0.0556016 -0.19630790 -5.3550160 -0.64571643
[23,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1916196 1.1447303 -3.087935 0.21698563 -0.0556016 4.0190707 1.63941245 1.8117645 1.34378117
[24,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.4950555 3.2432616 -1.486801 0.01800482 -0.1963079 1.6394124 2.75909588 -0.7756486 0.06479731
[25,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.6854842 -3.3624887 1.637401 1.32441392 -5.3550160 1.8117645 -0.77564863 10.0208796 3.29290470
[26,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.5653566 -0.7092751 1.303399 0.88865550 -0.6457164 1.3437812 0.06479731 3.2929047 1.31668449
> Hi = Ai + H0
> round(Hi, 3)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
[1,] 1.5 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.0 0 0 0.000 0.000 0.000 -1.000 0.000 0.000 0.000 0.000 0.000
[2,] 0.0 1.5 0.0 0.0 0.0 0.0 0.0 0.0 0 0.0 0.0 0.0 0.0 0.0 0.5 -1 0 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[3,] 0.5 0.0 1.5 0.0 0.0 0.0 0.0 0.0 0 0.0 0.0 0.0 0.0 0.0 0.0 0 0 0.000 0.000 0.000 -1.000 0.000 0.000 0.000 0.000 0.000
[4,] 0.0 0.0 0.0 1.5 0.0 0.0 0.0 0.0 0 0.0 0.0 0.0 0.5 0.0 -1.0 0 0 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[5,] 0.0 0.0 0.0 0.0 1.5 0.0 0.0 0.0 0 0.0 0.0 0.0 0.0 0.0 0.5 0 -1 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[6,] 0.0 0.0 0.0 0.0 0.0 1.5 0.0 0.0 0 0.0 0.0 0.0 0.0 0.5 0.0 0 0 -1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[7,] 0.0 0.0 0.0 0.0 0.0 0.0 1.5 0.0 0 0.0 0.0 0.0 0.0 0.5 0.0 0 0 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -1.000 0.000
[8,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.5 0 0.0 0.0 0.0 0.0 0.5 0.0 0 0 0.000 0.000 0.000 0.000 -1.000 0.000 0.000 0.000 0.000
[9,] 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 1.0 0.0 0 0 0.000 -1.000 -1.000 0.000 0.000 0.000 0.000 0.000 0.000
[10,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0 1.5 0.0 0.0 0.0 0.5 0.0 0 0 0.000 0.000 0.000 0.000 0.000 0.000 -1.000 0.000 0.000
[11,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0 0.0 1.5 0.0 0.0 0.5 0.0 0 0 0.000 0.000 0.000 0.000 0.000 -1.000 0.000 0.000 0.000
[12,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0 0.0 0.0 1.5 0.0 0.5 0.0 0 0 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -1.000
[13,] 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.0 0 0.0 0.0 0.0 1.5 0.0 -1.0 0 0 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[14,] 0.0 0.0 0.0 0.0 0.0 0.5 0.5 0.5 1 0.5 0.5 0.5 0.0 5.0 0.0 0 0 -1.000 -1.000 -1.000 0.000 -1.000 -1.000 -1.000 -1.000 -1.000
[15,] 0.0 0.5 0.0 -1.0 0.5 0.0 0.0 0.0 0 0.0 0.0 0.0 -1.0 0.0 3.0 -1 -1 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[16,] 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0 0.0 0.0 0.0 0.0 0.0 -1.0 2 0 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[17,] 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0 0.0 0.0 0.0 0.0 0.0 -1.0 0 2 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[18,] 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0 0.0 0.0 0.0 0.0 -1.0 0.0 0 0 3.329 1.029 1.002 0.410 1.649 0.192 1.495 -0.685 0.565
[19,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1 0.0 0.0 0.0 0.0 -1.0 0.0 0 0 1.029 6.206 -1.211 0.567 2.091 1.145 3.243 -3.362 -0.709
[20,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1 0.0 0.0 0.0 0.0 -1.0 0.0 0 0 1.002 -1.211 7.417 1.825 2.288 -3.088 -1.487 1.637 1.303
[21,] -1.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0 0.0 0.0 0.0 0.0 0.0 0.0 0 0 0.410 0.567 1.825 2.593 1.227 0.217 0.018 1.324 0.889
[22,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0 0.0 0.0 0.0 0.0 -1.0 0.0 0 0 1.649 2.091 2.288 1.227 8.289 -0.056 -0.196 -5.355 -0.646
[23,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0 0.0 -1.0 0.0 0.0 -1.0 0.0 0 0 0.192 1.145 -3.088 0.217 -0.056 6.019 1.639 1.812 1.344
[24,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0 -1.0 0.0 0.0 0.0 -1.0 0.0 0 0 1.495 3.243 -1.487 0.018 -0.196 1.639 4.759 -0.776 0.065
[25,] 0.0 0.0 0.0 0.0 0.0 0.0 -1.0 0.0 0 0.0 0.0 0.0 0.0 -1.0 0.0 0 0 -0.685 -3.362 1.637 1.324 -5.355 1.812 -0.776 12.021 3.293
[26,] 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0 0.0 0.0 -1.0 0.0 -1.0 0.0 0 0 0.565 -0.709 1.303 0.889 -0.646 1.344 0.065 3.293 3.317
>
> # weight R
> R = diag(data[1:10, 5])
> R
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 558 0 0 0 0 0 0 0 0 0
[2,] 0 722 0 0 0 0 0 0 0 0
[3,] 0 0 300 0 0 0 0 0 0 0
[4,] 0 0 0 73 0 0 0 0 0 0
[5,] 0 0 0 0 52 0 0 0 0 0
[6,] 0 0 0 0 0 87 0 0 0 0
[7,] 0 0 0 0 0 0 64 0 0 0
[8,] 0 0 0 0 0 0 0 103 0 0
[9,] 0 0 0 0 0 0 0 0 13 0
[10,] 0 0 0 0 0 0 0 0 0 125
> Ri = ginv(R)
> Ri
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0.001792115 0.000000000 0.000000000 0.00000000 0.00000000 0.00000000 0.000000 0.000000000 0.00000000 0.000
[2,] 0.000000000 0.001385042 0.000000000 0.00000000 0.00000000 0.00000000 0.000000 0.000000000 0.00000000 0.000
[3,] 0.000000000 0.000000000 0.003333333 0.00000000 0.00000000 0.00000000 0.000000 0.000000000 0.00000000 0.000
[4,] 0.000000000 0.000000000 0.000000000 0.01369863 0.00000000 0.00000000 0.000000 0.000000000 0.00000000 0.000
[5,] 0.000000000 0.000000000 0.000000000 0.00000000 0.01923077 0.00000000 0.000000 0.000000000 0.00000000 0.000
[6,] 0.000000000 0.000000000 0.000000000 0.00000000 0.00000000 0.01149425 0.000000 0.000000000 0.00000000 0.000
[7,] 0.000000000 0.000000000 0.000000000 0.00000000 0.00000000 0.00000000 0.015625 0.000000000 0.00000000 0.000
[8,] 0.000000000 0.000000000 0.000000000 0.00000000 0.00000000 0.00000000 0.000000 0.009708738 0.00000000 0.000
[9,] 0.000000000 0.000000000 0.000000000 0.00000000 0.00000000 0.00000000 0.000000 0.000000000 0.07692308 0.000
[10,] 0.000000000 0.000000000 0.000000000 0.00000000 0.00000000 0.00000000 0.000000 0.000000000 0.00000000 0.008
>
> # Left Hand Side
> LHS = rbind(
+ cbind(t(X) %*% Ri %*% X, t(X) %*% Ri %*% W),
+ cbind(t(W) %*% Ri %*% X, t(W) %*% Ri %*% W + Hi * alpha))
> round(LHS, 3)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26] [,27]
[1,] 0.161 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.002 0.001 0.003 0.014 0.019 0.011 0.016 0.010 0.077 0.008 0.000 0.000 0.000 0.000
[2,] 0.000 10.428 0.000 3.476 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.000 -6.952 0.000 0.000 0.000 0.000 0.000
[3,] 0.000 0.000 10.428 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 3.476 -6.952 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[4,] 0.000 3.476 0.000 10.428 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.000 -6.952 0.000 0.000 0.000 0.000 0.000
[5,] 0.000 0.000 0.000 0.000 10.428 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 3.476 0.000 -6.952 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[6,] 0.000 0.000 0.000 0.000 0.000 10.428 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 3.476 0.000 -6.952 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[7,] 0.000 0.000 0.000 0.000 0.000 0.000 10.428 0.000 0.000 0.000 0.000 0.000 0.000 0.000 3.476 0.000 0.000 0.000 -6.952 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[8,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 10.428 0.000 0.000 0.000 0.000 0.000 0.000 3.476 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -6.952 0.000
[9,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 10.428 0.000 0.000 0.000 0.000 0.000 3.476 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -6.952 0.000 0.000 0.000 0.000
[10,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 13.904 0.000 0.000 0.000 0.000 6.952 0.000 0.000 0.000 0.000 -6.952 -6.952 0.000 0.000 0.000 0.000 0.000 0.000
[11,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 10.428 0.000 0.000 0.000 3.476 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -6.952 0.000 0.000
[12,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 10.428 0.000 0.000 3.476 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -6.952 0.000 0.000 0.000
[13,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 10.428 0.000 3.476 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -6.952
[14,] 0.002 0.000 0.000 0.000 3.476 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 10.430 0.000 -6.952 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[15,] 0.001 0.000 0.000 0.000 0.000 0.000 3.476 3.476 3.476 6.952 3.476 3.476 3.476 0.000 34.762 0.000 0.000 0.000 -6.952 -6.952 -6.952 0.000 -6.952 -6.952 -6.952 -6.952 -6.952
[16,] 0.003 0.000 3.476 0.000 -6.952 3.476 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -6.952 0.000 20.860 -6.952 -6.952 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[17,] 0.014 0.000 -6.952 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -6.952 13.918 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[18,] 0.019 0.000 0.000 0.000 0.000 -6.952 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -6.952 0.000 13.923 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[19,] 0.011 0.000 0.000 0.000 0.000 0.000 -6.952 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -6.952 0.000 0.000 0.000 23.155 7.154 6.967 2.847 11.462 1.332 10.394 -4.766 3.930
[20,] 0.016 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -6.952 0.000 0.000 0.000 0.000 -6.952 0.000 0.000 0.000 7.154 43.157 -8.418 3.940 14.535 7.958 22.548 -23.376 -4.931
[21,] 0.010 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -6.952 0.000 0.000 0.000 0.000 -6.952 0.000 0.000 0.000 6.967 -8.418 51.576 12.685 15.910 -21.468 -10.336 11.383 9.061
[22,] 0.077 -6.952 0.000 -6.952 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 2.847 3.940 12.685 18.104 8.530 1.509 0.125 9.207 6.178
[23,] 0.008 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -6.952 0.000 0.000 0.000 0.000 0.000 -6.952 0.000 0.000 0.000 11.462 14.535 15.910 8.530 57.638 -0.387 -1.365 -37.229 -4.489
[24,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -6.952 0.000 0.000 -6.952 0.000 0.000 0.000 1.332 7.958 -21.468 1.509 -0.387 41.845 11.397 12.596 9.342
[25,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -6.952 0.000 0.000 0.000 -6.952 0.000 0.000 0.000 10.394 22.548 -10.336 0.125 -1.365 11.397 33.086 -5.392 0.450
[26,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -6.952 0.000 0.000 0.000 0.000 0.000 0.000 -6.952 0.000 0.000 0.000 -4.766 -23.376 11.383 9.207 -37.229 12.596 -5.392 83.571 22.893
[27,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -6.952 0.000 -6.952 0.000 0.000 0.000 3.930 -4.931 9.061 6.178 -4.489 9.342 0.450 22.893 23.058
>
> # Right Hand Side
> RHS = rbind(t(X) %*% Ri %*% y,
+ t(W) %*% Ri %*% y)
> round(RHS, 3)
[,1]
[1,] 1.351
[2,] 0.000
[3,] 0.000
[4,] 0.000
[5,] 0.000
[6,] 0.000
[7,] 0.000
[8,] 0.000
[9,] 0.000
[10,] 0.000
[11,] 0.000
[12,] 0.000
[13,] 0.000
[14,] 0.016
[15,] 0.019
[16,] 0.042
[17,] 0.211
[18,] 0.113
[19,] 0.089
[20,] 0.159
[21,] 0.047
[22,] 0.585
[23,] 0.070
[24,] 0.000
[25,] 0.000
[26,] 0.000
[27,] 0.000
>
> # Solutions
>
> # generalized inverse of LHS
> gi_LHS = ginv(LHS)
> round(gi_LHS, 3)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26] [,27]
[1,] 6.279 -0.080 -0.006 -0.080 -0.009 -0.009 -0.022 0.019 0.022 0.018 -0.006 -0.011 0.007 -0.011 0.026 -0.019 -0.019 -0.022 -0.020 0.036 0.025 -0.161 0.046 -0.004 0.004 0.042 0.023
[2,] -0.080 0.201 0.000 0.058 0.000 0.000 0.019 -0.022 -0.020 -0.029 0.016 0.008 -0.003 0.000 -0.031 0.000 0.000 0.000 0.013 -0.052 -0.036 0.187 -0.046 -0.003 0.008 -0.048 -0.020
[3,] -0.006 0.000 0.144 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.072 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[4,] -0.080 0.058 0.000 0.201 0.000 0.000 0.019 -0.022 -0.020 -0.029 0.016 0.008 -0.003 0.000 -0.031 0.000 0.000 0.000 0.013 -0.052 -0.036 0.187 -0.046 -0.003 0.008 -0.048 -0.020
[5,] -0.009 0.000 0.000 0.000 0.144 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.072 0.036 0.036 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[6,] -0.009 0.000 0.000 0.000 0.000 0.144 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.072 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[7,] -0.022 0.019 0.000 0.019 0.000 0.000 0.154 -0.004 -0.022 0.015 -0.026 0.022 -0.015 0.000 -0.021 0.000 0.000 0.000 0.077 0.007 0.003 0.038 -0.043 0.022 -0.050 -0.017 -0.033
[8,] 0.019 -0.022 0.000 -0.022 0.000 0.000 -0.004 0.128 0.031 -0.004 0.008 -0.020 0.000 0.000 -0.004 0.000 0.000 0.000 -0.009 0.010 -0.023 -0.043 0.045 -0.031 0.010 0.046 -0.002
[9,] 0.022 -0.020 0.000 -0.020 0.000 0.000 -0.022 0.031 0.150 -0.031 0.029 -0.037 0.021 0.000 -0.002 0.000 0.000 0.000 -0.034 -0.016 -0.048 -0.040 0.080 -0.057 0.042 0.046 0.030
[10,] 0.018 -0.029 0.000 -0.029 0.000 0.000 0.015 -0.004 -0.031 0.155 -0.037 0.050 -0.020 0.000 -0.015 0.000 0.000 0.000 0.015 0.068 0.084 -0.057 -0.054 0.067 -0.063 -0.014 -0.038
[11,] -0.006 0.016 0.000 0.016 0.000 0.000 -0.026 0.008 0.029 -0.037 0.161 -0.022 0.016 0.000 -0.015 0.000 0.000 0.000 -0.047 -0.056 -0.032 0.031 0.036 -0.041 0.090 0.005 0.017
[12,] -0.011 0.008 0.000 0.008 0.000 0.000 0.022 -0.020 -0.037 0.050 -0.022 0.162 -0.030 0.000 -0.020 0.000 0.000 0.000 0.023 0.013 0.066 0.016 -0.066 0.089 -0.043 -0.039 -0.056
[13,] 0.007 -0.003 0.000 -0.003 0.000 0.000 -0.015 0.000 0.021 -0.020 0.016 -0.030 0.164 0.000 -0.009 0.000 0.000 0.000 -0.028 -0.009 -0.041 -0.006 0.026 -0.050 0.020 -0.005 0.097
[14,] -0.011 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.144 0.000 0.072 0.036 0.036 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[15,] 0.026 -0.031 0.000 -0.031 0.000 0.000 -0.021 -0.004 -0.002 -0.015 -0.015 -0.020 -0.009 0.000 0.057 0.000 0.000 0.000 -0.003 0.017 0.009 -0.061 0.025 -0.001 0.006 0.022 0.014
[16,] -0.019 0.000 0.000 0.000 0.072 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.072 0.000 0.144 0.072 0.072 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[17,] -0.019 0.000 0.072 0.000 0.036 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.036 0.000 0.072 0.144 0.036 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[18,] -0.022 0.000 0.000 0.000 0.036 0.072 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.036 0.000 0.072 0.036 0.143 0.000 0.000 0.000 0.001 0.000 0.000 0.000 0.000 0.000
[19,] -0.020 0.013 0.000 0.013 0.000 0.000 0.077 -0.009 -0.034 0.015 -0.047 0.023 -0.028 0.000 -0.003 0.000 0.000 0.000 0.114 0.019 0.009 0.026 -0.052 0.033 -0.072 -0.014 -0.043
[20,] 0.036 -0.052 0.000 -0.052 0.000 0.000 0.007 0.010 -0.016 0.068 -0.056 0.013 -0.009 0.000 0.017 0.000 0.000 0.000 0.019 0.104 0.048 -0.104 -0.015 0.028 -0.076 0.023 -0.005
[21,] 0.025 -0.036 0.000 -0.036 0.000 0.000 0.003 -0.023 -0.048 0.084 -0.032 0.066 -0.041 0.000 0.009 0.000 0.000 0.000 0.009 0.048 0.128 -0.072 -0.067 0.104 -0.043 -0.029 -0.057
[22,] -0.161 0.187 0.000 0.187 0.000 0.000 0.038 -0.043 -0.040 -0.057 0.031 0.016 -0.006 0.000 -0.061 0.000 0.000 0.001 0.026 -0.104 -0.072 0.374 -0.091 -0.007 0.016 -0.095 -0.040
[23,] 0.046 -0.046 0.000 -0.046 0.000 0.000 -0.043 0.045 0.080 -0.054 0.036 -0.066 0.026 0.000 0.025 0.000 0.000 0.000 -0.052 -0.015 -0.067 -0.091 0.133 -0.086 0.066 0.080 0.052
[24,] -0.004 -0.003 0.000 -0.003 0.000 0.000 0.022 -0.031 -0.057 0.067 -0.041 0.089 -0.050 0.000 -0.001 0.000 0.000 0.000 0.033 0.028 0.104 -0.007 -0.086 0.133 -0.062 -0.048 -0.076
[25,] 0.004 0.008 0.000 0.008 0.000 0.000 -0.050 0.010 0.042 -0.063 0.090 -0.043 0.020 0.000 0.006 0.000 0.000 0.000 -0.072 -0.076 -0.043 0.016 0.066 -0.062 0.138 0.019 0.033
[26,] 0.042 -0.048 0.000 -0.048 0.000 0.000 -0.017 0.046 0.046 -0.014 0.005 -0.039 -0.005 0.000 0.022 0.000 0.000 0.000 -0.014 0.023 -0.029 -0.095 0.080 -0.048 0.019 0.081 0.004
[27,] 0.023 -0.020 0.000 -0.020 0.000 0.000 -0.033 -0.002 0.030 -0.038 0.017 -0.056 0.097 0.000 0.014 0.000 0.000 0.000 -0.043 -0.005 -0.057 -0.040 0.052 -0.076 0.033 0.004 0.152
>
> # solution
> sol = gi_LHS %*% RHS
>
> # general mean effect
> gme = sol[1]
> round(gme, 3)
[1] 8.39
>
> # (g)ebv of animals
> gebv = sol[2:27]
> round(gebv, 3)
[1] -0.012 0.007 -0.012 0.003 -0.003 -0.003 0.004 0.004 0.002 -0.002 -0.003 0.002 0.003 0.004 0.006 0.013 -0.002 -0.002 0.007 0.001 -0.024 0.008 -0.003 -0.001 0.008 0.005
>
결과가 책과 다르다. 책 187쪽에는 5세대의 혈통을 이용하여 구한 A(NRM)가 나오는데 책에서 제시한 혈통을 가지고는 이걸 구할 수가 없다. 본 예제에서 책 187쪽의 혈통을 사용한다고 하였는데 프로그램의 흐름상 사용할 수가 없었다. 이런 이유로 해서 같은 해가 나올 수 없는 것으로 보인다. 본 예제를 위의 참고 포스팅에서는 blupf90으로 푸는데 역시 해가 다르다. 그 이유는 참고 포스팅에서는 EDC를 가중치로 줄 때 그대로 사용하였다. EDC가 아니라 EDC의 역수를 가중치로 주어야 한다. EDC의 역수를 가중치로 주었을 때 동일한 해를 구할 수 있었다.
'Animal Breeding > R for Genetic Evaluation' 카테고리의 다른 글
| 혈통 그루핑 (0) | 2023.05.08 |
|---|---|
| R을 이용한 사회적 관계 모형(Social Interaction Model)의 육종가 구하기 (0) | 2021.04.15 |
| GBLUP Model with residual polygenic effects (0) | 2020.12.29 |
| SNP-BLUP Model with Polygenic Effects(unweighted analysis) (0) | 2020.12.24 |
| Mixed Linear Model(GBLUP model) for Computing DGV (0) | 2020.12.21 |