# GBLUP Model with residual polygenic effects
# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.
# Raphael Mrode
# Example 11.2 p189
GBLUP Model은 유전체 자료를 가지고 개체들의 DGV를 구하는 모형이다. 여기서 유전체 자료가 모든 상가적 유전 효과를 설명할 수 있다고 보기 어려워 polygenic effect를 추가한다. 간단한 설명은 다음 포스팅을 참고한다.
SNP-BLUP Model with Polygenic Effects(unweighted analysis)
# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition. # Raphael Mrode # Example 11.5 p189 SNP effect로 모든 상가적 유전 분산을 설명할 수 없을 수도 있다. 그래서 설명..
bhpark.tistory.com
data
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
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
program
# GBLUP Model with residual polygenic effects
# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.
# Raphael Mrode
# Example 11.2 p189
# 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) {
# 부모??? 모??? 경우
nrm[animal, animal] = 1
} else if (sire != 0 & dam == 0) {
# 부??? ?????? ?????? 경우
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) {
# 모만 ?????? ?????? 경우
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 {
# 부??? 모??? ?????? ?????? 경우
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("snp_data.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
data
# read pedigree : animal sire dam
pedi = read.table("snp_pedi.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
pedi
# variance component and ratio
sigma_a = 35.241
sigma_e = 245
sigma_a
sigma_e
# ratio of polygenic effect
rpe = 0.1
# variance of residual polygenic effect
sigma_u = sigma_a * rpe
sigma_u
# variance ratio, alpha
alpha1 = sigma_e / sigma_u
alpha1
alpha2 = sigma_e / (sigma_a - sigma_u)
alpha2
# SNP
M_all = as.matrix(data[, 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:8, 6]
y
# design matrix of fixed effect(general mean)
X = matrix(rep(1, 8), 8, 1)
X
# design matrix for polygenic effects
W1 = cbind(matrix(c(0), 8, 12), diag(8), matrix(c(0), 8, 6))
W1
# inverse matrix of NRM
Ai = ainv(pedi)
Ai
# design matrix for DGV
W2 = cbind(diag(8), matrix(c(0), 8, 6))
W2
# 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[13:26, 13:26]
A22
# weighted G
wt = 0.99
Gw = wt * G + (1 - wt) * A22
round(Gw, 3)
Gwi = ginv(Gw)
round(Gwi, 3)
# Left Hand Side
LHS = rbind(
cbind(t(X) %*% X, t(X) %*% W1, t(X) %*% W2),
cbind(t(W1) %*% X, t(W1) %*% W1 + Ai * alpha1, t(W1) %*% W2),
cbind(t(W2) %*% X, t(W2) %*% W1, t(W2) %*% W2 + Gwi * alpha2))
round(LHS, 3)
# Right Hand Side
RHS = rbind(t(X) %*% y,
t(W1) %*% y,
t(W2) %*% y)
round(RHS, 3)
# Solutions
# generalized inverse of LHS
gi_LHS = ginv(LHS)
round(gi_LHS, 3)
# solution
sol = gi_LHS %*% RHS
round(sol, 3)
RStudio에서 실행한 화면
프로그램 실행 결과
> # GBLUP Model with residual polygenic effects
>
> # Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.
>
> # Raphael Mrode
>
> # Example 11.2 p189
>
> # 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) {
+ # 부모??? 모??? 경우
+ nrm[animal, animal] = 1
+ } else if (sire != 0 & dam == 0) {
+ # 부??? ?????? ?????? 경우
+ 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) {
+ # 모만 ?????? ?????? 경우
+ 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 {
+ # 부??? 모??? ?????? ?????? 경우
+ 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/2020/job/공부하기/07_Linear Models for the Prediction of Animal Breeding Values_3rd_Mrode/23_GBLUP with polygenic effect"
>
> # read data
> data = read.table("snp_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("snp_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
>
> # ratio of polygenic effect
> rpe = 0.1
>
> # variance of residual polygenic effect
> sigma_u = sigma_a * rpe
> sigma_u
[1] 3.5241
>
> # variance ratio, alpha
> alpha1 = sigma_e / sigma_u
> alpha1
[1] 69.5213
>
> alpha2 = sigma_e / (sigma_a - sigma_u)
> alpha2
[1] 7.724588
>
> # SNP
> M_all = as.matrix(data[, 7:16])
> M_all
snp1 snp2 snp3 snp4 snp5 snp6 snp7 snp8 snp9 snp10
[1,] 2 0 1 1 0 0 0 2 1 2
[2,] 1 0 0 0 0 2 0 2 1 0
[3,] 1 1 2 1 1 0 0 2 1 2
[4,] 0 0 2 1 0 1 0 2 2 1
[5,] 0 1 1 2 0 0 0 2 1 2
[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.32142857 0.17857143 0.35714286 0.35714286 0.14285714 0.60714286 0.07142857 0.96428571 0.57142857 0.39285714
>
> # sum of 2pq(heterozygote ratio)
> sum2pq = sum(2 * snp_mean * (1 - snp_mean))
> sum2pq
[1] 3.538265
>
> # observation
> y = data[1:8, 6]
> y
[1] 9.0 13.4 12.7 15.4 5.9 7.7 10.2 4.8
>
> # design matrix of fixed effect(general mean)
> X = matrix(rep(1, 8), 8, 1)
> X
[,1]
[1,] 1
[2,] 1
[3,] 1
[4,] 1
[5,] 1
[6,] 1
[7,] 1
[8,] 1
>
> # design matrix for polygenic effects
> W1 = cbind(matrix(c(0), 8, 12), diag(8), matrix(c(0), 8, 6))
> W1
[,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
> # 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
>
> # design matrix for DGV
> W2 = cbind(diag(8), matrix(c(0), 8, 6))
> W2
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 1 0 0 0 0 0 0 0 0 0 0 0 0 0
[2,] 0 1 0 0 0 0 0 0 0 0 0 0 0 0
[3,] 0 0 1 0 0 0 0 0 0 0 0 0 0 0
[4,] 0 0 0 1 0 0 0 0 0 0 0 0 0 0
[5,] 0 0 0 0 1 0 0 0 0 0 0 0 0 0
[6,] 0 0 0 0 0 1 0 0 0 0 0 0 0 0
[7,] 0 0 0 0 0 0 1 0 0 0 0 0 0 0
[8,] 0 0 0 0 0 0 0 1 0 0 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
[1,] 1.357 -0.357 0.286 0.286 -0.286 -1.214 -0.143 0.071 -0.143 1.214
[2,] 0.357 -0.357 -0.714 -0.714 -0.286 0.786 -0.143 0.071 -0.143 -0.786
[3,] 0.357 0.643 1.286 0.286 0.714 -1.214 -0.143 0.071 -0.143 1.214
[4,] -0.643 -0.357 1.286 0.286 -0.286 -0.214 -0.143 0.071 0.857 0.214
[5,] -0.643 0.643 0.286 1.286 -0.286 -1.214 -0.143 0.071 -0.143 1.214
[6,] 0.357 0.643 -0.714 0.286 -0.286 0.786 -0.143 0.071 0.857 0.214
[7,] -0.643 -0.357 0.286 0.286 -0.286 0.786 -0.143 0.071 0.857 -0.786
[8,] -0.643 0.643 0.286 -0.714 -0.286 -0.214 -0.143 0.071 0.857 -0.786
[9,] 1.357 -0.357 -0.714 -0.714 -0.286 -0.214 1.857 0.071 -0.143 1.214
[10,] -0.643 -0.357 -0.714 0.286 0.714 0.786 -0.143 0.071 -1.143 -0.786
[11,] -0.643 0.643 0.286 -0.714 -0.286 -0.214 -0.143 0.071 0.857 0.214
[12,] 0.357 -0.357 -0.714 -0.714 0.714 -0.214 -0.143 0.071 -1.143 -0.786
[13,] -0.643 -0.357 -0.714 0.286 0.714 0.786 -0.143 0.071 -0.143 -0.786
[14,] 0.357 -0.357 0.286 0.286 -0.286 0.786 -0.143 -0.929 -1.143 -0.786
> G = Z %*% t(Z) / sum2pq
> round(G, 3)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 1.472 -0.446 0.988 0.059 0.685 -0.163 -0.708 -0.547 0.887 -0.789 -0.203 -0.143 -0.829 -0.264
[2,] -0.446 0.745 -0.930 -0.446 -0.950 0.180 0.200 0.079 0.099 0.402 -0.143 0.483 0.362 0.362
[3,] 0.988 -0.930 1.634 0.422 1.048 -0.365 -0.627 -0.183 0.120 -0.708 0.160 -0.345 -0.748 -0.466
[4,] 0.059 -0.446 0.422 0.907 0.402 -0.163 0.422 0.301 -0.526 -0.506 0.362 -0.708 -0.264 -0.264
[5,] 0.685 -0.950 1.048 0.402 1.593 -0.102 -0.365 -0.203 -0.183 -0.446 0.140 -0.647 -0.486 -0.486
[6,] -0.163 0.180 -0.365 -0.163 -0.102 0.745 0.200 0.079 0.099 -0.163 0.140 -0.365 0.079 -0.203
[7,] -0.708 0.200 -0.627 0.422 -0.365 0.200 0.786 0.382 -0.728 0.140 0.160 -0.345 0.382 0.099
[8,] -0.547 0.079 -0.183 0.301 -0.203 0.079 0.382 0.826 -0.567 -0.264 0.604 -0.183 -0.022 -0.304
[9,] 0.887 0.099 0.120 -0.526 -0.183 0.099 -0.728 -0.567 2.280 -0.526 -0.224 0.120 -0.567 -0.284
[10,] -0.789 0.402 -0.708 -0.506 -0.446 -0.163 0.140 -0.264 -0.526 1.190 -0.486 0.705 0.867 0.584
[11,] -0.203 -0.143 0.160 0.362 0.140 0.140 0.160 0.604 -0.224 -0.486 0.665 -0.405 -0.244 -0.526
[12,] -0.143 0.483 -0.345 -0.708 -0.647 -0.365 -0.345 -0.183 0.120 0.705 -0.405 1.068 0.382 0.382
[13,] -0.829 0.362 -0.748 -0.264 -0.486 0.079 0.382 -0.022 -0.567 0.867 -0.244 0.382 0.826 0.261
[14,] -0.264 0.362 -0.466 -0.264 -0.486 -0.203 0.099 -0.304 -0.284 0.584 -0.526 0.382 0.261 1.109
>
> # 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[13:26, 13:26]
> A22
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 1.00 0.0 0.5 0.25 0.25 0.00 0.00 0.00 0 0.00 0.00 0.00 0.00 0.00
[2,] 0.00 1.0 0.0 0.00 0.00 0.50 0.50 0.50 0 0.50 0.50 0.50 0.50 0.50
[3,] 0.50 0.0 1.0 0.50 0.50 0.00 0.00 0.00 0 0.00 0.00 0.00 0.00 0.00
[4,] 0.25 0.0 0.5 1.00 0.25 0.00 0.00 0.00 0 0.00 0.00 0.00 0.00 0.00
[5,] 0.25 0.0 0.5 0.25 1.00 0.00 0.00 0.00 0 0.00 0.00 0.00 0.00 0.00
[6,] 0.00 0.5 0.0 0.00 0.00 1.00 0.25 0.25 0 0.25 0.25 0.25 0.25 0.25
[7,] 0.00 0.5 0.0 0.00 0.00 0.25 1.00 0.50 0 0.25 0.25 0.25 0.25 0.25
[8,] 0.00 0.5 0.0 0.00 0.00 0.25 0.50 1.00 0 0.25 0.25 0.25 0.25 0.25
[9,] 0.00 0.0 0.0 0.00 0.00 0.00 0.00 0.00 1 0.00 0.00 0.00 0.00 0.00
[10,] 0.00 0.5 0.0 0.00 0.00 0.25 0.25 0.25 0 1.00 0.25 0.25 0.25 0.25
[11,] 0.00 0.5 0.0 0.00 0.00 0.25 0.25 0.25 0 0.25 1.00 0.25 0.25 0.25
[12,] 0.00 0.5 0.0 0.00 0.00 0.25 0.25 0.25 0 0.25 0.25 1.00 0.25 0.25
[13,] 0.00 0.5 0.0 0.00 0.00 0.25 0.25 0.25 0 0.25 0.25 0.25 1.00 0.25
[14,] 0.00 0.5 0.0 0.00 0.00 0.25 0.25 0.25 0 0.25 0.25 0.25 0.25 1.00
>
> # weighted G
> wt = 0.99
> Gw = wt * G + (1 - wt) * A22
> round(Gw, 3)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 1.468 -0.441 0.983 0.061 0.681 -0.161 -0.701 -0.541 0.878 -0.781 -0.201 -0.141 -0.821 -0.261
[2,] -0.441 0.748 -0.921 -0.441 -0.941 0.183 0.203 0.084 0.099 0.403 -0.136 0.483 0.363 0.363
[3,] 0.983 -0.921 1.627 0.423 1.043 -0.361 -0.621 -0.181 0.118 -0.701 0.158 -0.341 -0.741 -0.461
[4,] 0.061 -0.441 0.423 0.908 0.401 -0.161 0.418 0.298 -0.521 -0.501 0.358 -0.701 -0.261 -0.261
[5,] 0.681 -0.941 1.043 0.401 1.587 -0.101 -0.361 -0.201 -0.181 -0.441 0.138 -0.641 -0.481 -0.481
[6,] -0.161 0.183 -0.361 -0.161 -0.101 0.748 0.201 0.081 0.099 -0.159 0.141 -0.359 0.081 -0.199
[7,] -0.701 0.203 -0.621 0.418 -0.361 0.201 0.788 0.383 -0.721 0.141 0.161 -0.339 0.381 0.101
[8,] -0.541 0.084 -0.181 0.298 -0.201 0.081 0.383 0.828 -0.561 -0.259 0.601 -0.179 -0.019 -0.299
[9,] 0.878 0.099 0.118 -0.521 -0.181 0.099 -0.721 -0.561 2.267 -0.521 -0.221 0.118 -0.561 -0.281
[10,] -0.781 0.403 -0.701 -0.501 -0.441 -0.159 0.141 -0.259 -0.521 1.188 -0.479 0.701 0.860 0.581
[11,] -0.201 -0.136 0.158 0.358 0.138 0.141 0.161 0.601 -0.221 -0.479 0.668 -0.399 -0.239 -0.519
[12,] -0.141 0.483 -0.341 -0.701 -0.641 -0.359 -0.339 -0.179 0.118 0.701 -0.399 1.068 0.381 0.381
[13,] -0.821 0.363 -0.741 -0.261 -0.481 0.081 0.381 -0.019 -0.561 0.860 -0.239 0.381 0.828 0.261
[14,] -0.261 0.363 -0.461 -0.261 -0.481 -0.199 0.101 -0.299 -0.281 0.581 -0.519 0.381 0.261 1.108
> Gwi = ginv(Gw)
> round(Gwi, 3)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 36.374 -27.981 -8.256 -7.273 -16.994 1.865 1.034 23.580 2.958 39.443 15.760 -17.655 0.710 5.177
[2,] -27.981 88.271 31.300 -7.935 26.742 -30.838 -16.075 -3.254 -1.422 -44.402 -27.862 -31.454 38.437 -6.093
[3,] -8.256 31.300 17.521 -13.840 10.826 -13.961 10.967 -4.374 1.150 -15.361 -3.670 -9.610 12.371 -1.151
[4,] -7.273 -7.935 -13.840 59.030 1.183 35.036 -44.042 13.642 4.069 8.965 -16.849 24.984 -3.933 4.327
[5,] -16.994 26.742 10.826 1.183 16.560 -3.741 6.055 -10.356 1.661 -23.795 -5.410 8.979 8.462 0.465
[6,] 1.865 -30.838 -13.961 35.036 -3.741 33.909 -9.403 4.611 4.196 17.932 1.618 33.086 -19.287 5.158
[7,] 1.034 -16.075 10.967 -44.042 6.055 -9.403 76.599 -30.644 2.333 -4.681 31.963 24.525 -23.251 0.443
[8,] 23.580 -3.254 -4.374 13.642 -10.356 4.611 -30.644 38.819 3.866 31.071 -11.692 -26.201 11.746 3.079
[9,] 2.958 -1.422 1.150 4.069 1.661 4.196 2.333 3.866 3.572 4.523 2.296 3.933 1.097 2.907
[10,] 39.443 -44.402 -15.361 8.965 -23.795 17.932 -4.681 31.071 4.523 65.241 12.420 -9.861 -22.528 2.734
[11,] 15.760 -27.862 -3.670 -16.849 -5.410 1.618 31.963 -11.692 2.296 12.420 33.290 12.724 -8.689 6.556
[12,] -17.655 -31.454 -9.610 24.984 8.979 33.086 24.525 -26.201 3.933 -9.861 12.724 65.530 -29.485 5.120
[13,] 0.710 38.437 12.371 -3.933 8.462 -19.287 -23.251 11.746 1.097 -22.528 -8.689 -29.485 46.431 4.505
[14,] 5.177 -6.093 -1.151 4.327 0.465 5.158 0.443 3.079 2.907 2.734 6.556 5.120 4.505 6.133
>
> # Left Hand Side
> LHS = rbind(
+ cbind(t(X) %*% X, t(X) %*% W1, t(X) %*% W2),
+ cbind(t(W1) %*% X, t(W1) %*% W1 + Ai * alpha1, t(W1) %*% W2),
+ cbind(t(W2) %*% X, t(W2) %*% W1, t(W2) %*% W2 + Gwi * alpha2))
> 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]
[1,] 8 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.000 0.000
[2,] 0 104.282 0.000 34.761 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 -69.521 0.000
[3,] 0 0.000 104.282 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 34.761 -69.521 0.000 0.000 0.000 0.000 0.000 0.000
[4,] 0 34.761 0.000 104.282 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 -69.521 0.000
[5,] 0 0.000 0.000 0.000 104.282 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 34.761 0.000 -69.521 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[6,] 0 0.000 0.000 0.000 0.000 104.282 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 34.761 0.000 -69.521 0.000 0.000 0.000 0.000 0.000
[7,] 0 0.000 0.000 0.000 0.000 0.000 104.282 0.000 0.000 0.000 0.000 0.000 0.000 0.000 34.761 0.000 0.000 0.000 -69.521 0.000 0.000 0.000 0.000
[8,] 0 0.000 0.000 0.000 0.000 0.000 0.000 104.282 0.000 0.000 0.000 0.000 0.000 0.000 34.761 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[9,] 0 0.000 0.000 0.000 0.000 0.000 0.000 0.000 104.282 0.000 0.000 0.000 0.000 0.000 34.761 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -69.521
[10,] 0 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 139.043 0.000 0.000 0.000 0.000 69.521 0.000 0.000 0.000 0.000 -69.521 -69.521 0.000 0.000
[11,] 0 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 104.282 0.000 0.000 0.000 34.761 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[12,] 0 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 104.282 0.000 0.000 34.761 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[13,] 0 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 104.282 0.000 34.761 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[14,] 1 0.000 0.000 0.000 34.761 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 105.282 0.000 -69.521 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[15,] 1 0.000 0.000 0.000 0.000 0.000 34.761 34.761 34.761 69.521 34.761 34.761 34.761 0.000 348.606 0.000 0.000 0.000 -69.521 -69.521 -69.521 0.000 -69.521
[16,] 1 0.000 34.761 0.000 -69.521 34.761 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -69.521 0.000 209.564 -69.521 -69.521 0.000 0.000 0.000 0.000 0.000
[17,] 1 0.000 -69.521 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -69.521 140.043 0.000 0.000 0.000 0.000 0.000 0.000
[18,] 1 0.000 0.000 0.000 0.000 -69.521 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -69.521 0.000 140.043 0.000 0.000 0.000 0.000 0.000
[19,] 1 0.000 0.000 0.000 0.000 0.000 -69.521 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -69.521 0.000 0.000 0.000 140.043 0.000 0.000 0.000 0.000
[20,] 1 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -69.521 0.000 0.000 0.000 0.000 -69.521 0.000 0.000 0.000 0.000 140.043 0.000 0.000 0.000
[21,] 1 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -69.521 0.000 0.000 0.000 0.000 -69.521 0.000 0.000 0.000 0.000 0.000 140.043 0.000 0.000
[22,] 0 -69.521 0.000 -69.521 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 139.043 0.000
[23,] 0 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -69.521 0.000 0.000 0.000 0.000 0.000 -69.521 0.000 0.000 0.000 0.000 0.000 0.000 0.000 139.043
[24,] 0 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 -69.521 0.000 0.000 -69.521 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[,24] [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41]
[1,] 0.000 0.000 0.000 0.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000
[2,] 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 0.000
[3,] 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 0.000
[4,] 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 0.000
[5,] 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 0.000
[6,] 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 0.000
[7,] 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 0.000
[8,] 0.000 0.000 -69.521 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
[9,] 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 0.000
[10,] 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 0.000
[11,] 0.000 -69.521 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
[12,] -69.521 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
[13,] 0.000 0.000 0.000 -69.521 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
[14,] 0.000 0.000 0.000 0.000 1.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
[15,] -69.521 -69.521 -69.521 -69.521 0.000 1.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
[16,] 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[17,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[18,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[19,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[20,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[21,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.000 0.000 0.000 0.000 0.000 0.000
[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.000 0.000
[23,] 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 0.000
[24,] 139.043 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
[ getOption("max.print") 에 도달했습니다 -- 17 행들을 생략합니다 ]
>
> # Right Hand Side
> RHS = rbind(t(X) %*% y,
+ t(W1) %*% y,
+ t(W2) %*% y)
> round(RHS, 3)
[,1]
[1,] 79.1
[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,] 9.0
[15,] 13.4
[16,] 12.7
[17,] 15.4
[18,] 5.9
[19,] 7.7
[20,] 10.2
[21,] 4.8
[22,] 0.0
[23,] 0.0
[24,] 0.0
[25,] 0.0
[26,] 0.0
[27,] 0.0
[28,] 9.0
[29,] 13.4
[30,] 12.7
[31,] 15.4
[32,] 5.9
[33,] 7.7
[34,] 10.2
[35,] 4.8
[36,] 0.0
[37,] 0.0
[38,] 0.0
[39,] 0.0
[40,] 0.0
[41,] 0.0
>
> # 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,] 0.142 0.000 -0.001 0.000 -0.002 -0.001 -0.001 0.000 0.000 -0.002 0.000 0.000 0.000 -0.003 -0.005 -0.004 -0.003 -0.003 -0.004 -0.004 -0.004 0.000 -0.002 -0.002 -0.002 -0.002 -0.002
[2,] 0.000 0.014 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 0.000 0.000 0.007 0.000 0.000 0.000 0.000 0.000
[3,] -0.001 0.000 0.014 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.007 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[4,] 0.000 0.000 0.000 0.014 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 0.007 0.000 0.000 0.000 0.000 0.000
[5,] -0.002 0.000 0.000 0.000 0.014 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.007 0.004 0.004 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[6,] -0.001 0.000 0.000 0.000 0.000 0.014 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.007 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[7,] -0.001 0.000 0.000 0.000 0.000 0.000 0.014 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.007 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 0.014 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 0.007 0.000
[9,] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.014 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.007 0.000 0.000 0.000 0.000
[10,] -0.002 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.014 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.007 0.007 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 0.014 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.007 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 0.014 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.007 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 0.014 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.007
[14,] -0.003 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.014 0.000 0.007 0.004 0.004 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[15,] -0.005 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.014 0.000 0.000 0.000 0.007 0.007 0.007 0.000 0.007 0.007 0.007 0.007 0.007
[16,] -0.004 0.000 0.000 0.000 0.007 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.007 0.000 0.014 0.007 0.007 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[17,] -0.003 0.000 0.007 0.000 0.004 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.004 0.000 0.007 0.014 0.004 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[18,] -0.003 0.000 0.000 0.000 0.004 0.007 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.004 0.000 0.007 0.004 0.014 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[19,] -0.004 0.000 0.000 0.000 0.000 0.000 0.007 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.007 0.000 0.000 0.000 0.014 0.004 0.004 0.000 0.004 0.004 0.004 0.004 0.004
[20,] -0.004 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.007 0.000 0.000 0.000 0.000 0.007 0.000 0.000 0.000 0.004 0.014 0.007 0.000 0.004 0.004 0.004 0.004 0.004
[21,] -0.004 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.007 0.000 0.000 0.000 0.000 0.007 0.000 0.000 0.000 0.004 0.007 0.014 0.000 0.004 0.004 0.004 0.004 0.004
[22,] 0.000 0.007 0.000 0.007 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 0.014 0.000 0.000 0.000 0.000 0.000
[23,] -0.002 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.007 0.000 0.000 0.000 0.000 0.000 0.007 0.000 0.000 0.000 0.004 0.004 0.004 0.000 0.014 0.004 0.004 0.004 0.004
[24,] -0.002 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.007 0.000 0.000 0.007 0.000 0.000 0.000 0.004 0.004 0.004 0.000 0.004 0.014 0.004 0.004 0.004
[,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36] [,37] [,38] [,39] [,40] [,41]
[1,] -0.017 0.018 -0.023 -0.027 -0.025 -0.009 -0.007 -0.013 0.012 0.033 -0.018 0.032 0.020 0.021
[2,] 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
[3,] 0.000 0.000 0.000 -0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[4,] 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
[5,] -0.001 0.001 -0.001 0.000 -0.001 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.001 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.001 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 0.000 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.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
[10,] 0.001 0.000 0.001 0.000 0.001 0.000 -0.001 -0.001 0.001 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 0.000 0.000 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 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 0.000 0.000
[14,] -0.002 0.001 -0.002 0.000 -0.001 0.001 0.001 0.001 -0.001 0.001 0.000 0.000 0.001 0.000
[15,] 0.002 -0.002 0.003 0.001 0.002 -0.001 -0.001 -0.001 0.000 -0.001 0.000 -0.001 -0.001 -0.001
[16,] -0.002 0.002 -0.003 -0.001 -0.002 0.001 0.001 0.001 0.000 0.001 0.000 0.000 0.001 0.001
[17,] -0.001 0.001 -0.002 -0.001 -0.001 0.001 0.000 0.000 0.000 0.001 0.000 0.001 0.001 0.000
[18,] -0.001 0.002 -0.002 0.000 -0.002 0.001 0.001 0.001 0.000 0.001 0.000 0.001 0.001 0.001
[19,] 0.001 -0.001 0.002 0.001 0.001 -0.001 -0.001 0.000 0.000 -0.001 0.000 0.000 -0.001 0.000
[20,] 0.002 -0.001 0.002 0.000 0.002 0.000 -0.001 -0.001 0.001 -0.001 0.000 0.000 -0.001 -0.001
[21,] 0.002 -0.001 0.002 0.000 0.002 0.000 -0.001 -0.001 0.001 -0.001 0.000 -0.001 -0.001 0.000
[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
[23,] 0.001 -0.001 0.001 0.001 0.001 0.000 0.000 0.000 0.000 -0.001 0.000 -0.001 -0.001 0.000
[24,] 0.001 -0.001 0.001 0.001 0.001 0.000 0.000 0.000 0.000 -0.001 0.000 -0.001 -0.001 0.000
[ getOption("max.print") 에 도달했습니다 -- 17 행들을 생략합니다 ]
>
> # solution
> sol = gi_LHS %*% RHS
> round(sol, 3)
[,1]
[1,] 9.936
[2,] 0.000
[3,] 0.037
[4,] 0.000
[5,] 0.025
[6,] -0.026
[7,] -0.014
[8,] 0.000
[9,] 0.000
[10,] -0.034
[11,] 0.000
[12,] 0.000
[13,] 0.000
[14,] 0.011
[15,] 0.001
[16,] 0.043
[17,] 0.077
[18,] -0.017
[19,] -0.020
[20,] -0.016
[21,] -0.052
[22,] 0.000
[23,] 0.000
[24,] 0.000
[25,] 0.000
[26,] 0.000
[27,] 0.000
[28,] 0.055
[29,] 0.107
[30,] 0.032
[31,] 0.228
[32,] -0.454
[33,] -0.317
[34,] 0.132
[35,] -0.201
[36,] 0.023
[37,] 0.107
[38,] -0.214
[39,] 0.132
[40,] 0.053
[41,] 0.318
>
'Animal Breeding > R for Genetic Evaluation' 카테고리의 다른 글
| R을 이용한 사회적 관계 모형(Social Interaction Model)의 육종가 구하기 (0) | 2021.04.15 |
|---|---|
| Single-step GBLUP(ssGBLUP) (0) | 2021.01.05 |
| SNP-BLUP Model with Polygenic Effects(unweighted analysis) (0) | 2020.12.24 |
| Mixed Linear Model(GBLUP model) for Computing DGV (0) | 2020.12.21 |
| Mixed Linear Model(SNP-BLUP Model) for Computing SNP effect, weighted analysis (0) | 2020.12.19 |