# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.
# Raphael Mrode
# Example 11.1 p180
모형에 대한 간단한 설명은 아래 글을 참조한다.
Fixed Effect Model for SNP Effects, unweighted analysis
# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition. # Raphael Mrode # Example 11.1 p180 SNP effect를 고정효과로 모형에 넣어 SNP 추정하고 개체의 genotype을 이용..
bhpark.tistory.com
여기서는 수소(bull)의 EDC(effective daughter contribution)를 가중치로 주어 모형을 분석한다. 책의 MME에 약간의 오타가 나오는데 아래 프로그램을 보고 올바른 MME를 유추하길 바란다.
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
# Fixed Effect Model for SNP Effect, weighted analysis
# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.
# Raphael Mrode
# Example 11.1 p180
# MASS packages
if (!("MASS" %in% installed.packages())) {
install.packages("MASS", dependencies = TRUE)
}
library(MASS)
# 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(".")
# print working directory
getwd()
# prdigree and data
# read pedigree : animal sire dam
pedi = read.table("fem_snp_pedi.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
# print
pedi
# read data : animal, sex, pre-weaning gain
data = read.table("fem_snp_data.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
# print
data
# variance component and ratio
sigma_a = 35.241
sigma_e = 245
alpha = sigma_e/sigma_a
# print
sigma_a
sigma_e
alpha
# observation
y = data[1:8, 6]
# print
y
# design matrix
# design matrix of fixed effect
X = matrix(rep(1, 8), 8, 1)
# print
X
# SNP effect
M = as.matrix(data[1:8, 7:9])
# design matrix for SNP effect
# mean of SNP effect
M_all = as.matrix(data[, 7:9])
snp_mean = colMeans(M_all) / 2
Z = sweep(M, 2, snp_mean * 2)
Z
# design matrix for polygenic effect
W = cbind(matrix(rep(0, 8 * 12), 8, 12), diag(8), matrix(rep(0, 8 * 6), 8, 6))
W
# Weight : R from EDC
R = diag(data[1:8, 5])
R
Rinv = ginv(R)
Rinv
# inverse matrix of NRM
ai = ainv(pedi)
# print
ai
# LHS, RHS
# LHS
LHS = rbind(
cbind(t(X) %*% Rinv %*% X, t(X) %*% Rinv %*% Z, t(X) %*% Rinv %*% W),
cbind(t(Z) %*% Rinv %*% X, t(Z) %*% Rinv %*% Z, t(Z) %*% Rinv %*% W),
cbind(t(W) %*% Rinv %*% X, t(W) %*% Rinv %*% Z, t(W) %*% Rinv %*% W + ai * alpha))
# print
LHS
# RHS
RHS = rbind(t(X) %*% Rinv %*% y,
t(Z) %*% Rinv %*% y,
t(W) %*% Rinv %*% y)
# print
RHS
# Solutions
# generalized inverse of LHS
gi_LHS = ginv(LHS)
# print
gi_LHS
# solution
sol = gi_LHS %*% RHS
# print
sol
# mean effect of solutions
mean_eff = sol[1]
mean_eff
# snp effect of solutions
snp_eff = sol[2:4]
snp_eff
# polygenic effect of solutions
polge_eff = sol[5:30]
polge_eff
# calculate DGV
# genotype of animals
# design matrix of animals
P2 = sweep(M_all, 2, snp_mean * 2)
P2
# DGV
dgv = P2 %*% snp_eff
dgv
# GEBV : animal 13 ~ 26
gebv = dgv + polge_eff[13:26]
gebv
실행
Result
> # Fixed Effect Model for SNP Effect, weighted analysis
>
> # Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.
>
> # Raphael Mrode
>
> # Example 11.1 p180
>
> # MASS packages
> if (!("MASS" %in% installed.packages())) {
+ install.packages("MASS", dependencies = TRUE)
+ }
> library(MASS)
>
> # 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(".")
>
> # print working directory
> getwd()
[1] "D:/users/bhpark/2020/job/공부하기/07_Linear Models for the Prediction of Animal Breeding Values_3rd_Mrode/12_fixed effecct model for SNP effects_weighted"
>
> # prdigree and data
>
> # read pedigree : animal sire dam
> pedi = read.table("fem_snp_pedi.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
>
> # print
> 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
>
> # read data : animal, sex, pre-weaning gain
> data = read.table("fem_snp_data.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
>
> # print
> 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
>
> # variance component and ratio
> sigma_a = 35.241
> sigma_e = 245
> alpha = sigma_e/sigma_a
>
> # print
> sigma_a
[1] 35.241
> sigma_e
[1] 245
> alpha
[1] 6.95213
>
> # observation
> y = data[1:8, 6]
>
> # print
> y
[1] 9.0 13.4 12.7 15.4 5.9 7.7 10.2 4.8
>
> # design matrix
>
> # design matrix of fixed effect
> X = matrix(rep(1, 8), 8, 1)
>
> # print
> X
[,1]
[1,] 1
[2,] 1
[3,] 1
[4,] 1
[5,] 1
[6,] 1
[7,] 1
[8,] 1
>
> # SNP effect
> M = as.matrix(data[1:8, 7:9])
>
> # design matrix for SNP effect
> # mean of SNP effect
> M_all = as.matrix(data[, 7:9])
> snp_mean = colMeans(M_all) / 2
> Z = sweep(M, 2, snp_mean * 2)
> Z
snp1 snp2 snp3
1 1.3571429 -0.3571429 0.2857143
2 0.3571429 -0.3571429 -0.7142857
3 0.3571429 0.6428571 1.2857143
4 -0.6428571 -0.3571429 1.2857143
5 -0.6428571 0.6428571 0.2857143
6 0.3571429 0.6428571 -0.7142857
7 -0.6428571 -0.3571429 0.2857143
8 -0.6428571 0.6428571 0.2857143
>
> # design matrix for polygenic effect
> W = cbind(matrix(rep(0, 8 * 12), 8, 12), diag(8), matrix(rep(0, 8 * 6), 8, 6))
> 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
>
> # Weight : R from EDC
> R = diag(data[1:8, 5])
> R
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 558 0 0 0 0 0 0 0
[2,] 0 722 0 0 0 0 0 0
[3,] 0 0 300 0 0 0 0 0
[4,] 0 0 0 73 0 0 0 0
[5,] 0 0 0 0 52 0 0 0
[6,] 0 0 0 0 0 87 0 0
[7,] 0 0 0 0 0 0 64 0
[8,] 0 0 0 0 0 0 0 103
> Rinv = ginv(R)
> Rinv
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 0.001792115 0.000000000 0.000000000 0.00000000 0.00000000 0.00000000 0.000000 0.000000000
[2,] 0.000000000 0.001385042 0.000000000 0.00000000 0.00000000 0.00000000 0.000000 0.000000000
[3,] 0.000000000 0.000000000 0.003333333 0.00000000 0.00000000 0.00000000 0.000000 0.000000000
[4,] 0.000000000 0.000000000 0.000000000 0.01369863 0.00000000 0.00000000 0.000000 0.000000000
[5,] 0.000000000 0.000000000 0.000000000 0.00000000 0.01923077 0.00000000 0.000000 0.000000000
[6,] 0.000000000 0.000000000 0.000000000 0.00000000 0.00000000 0.01149425 0.000000 0.000000000
[7,] 0.000000000 0.000000000 0.000000000 0.00000000 0.00000000 0.00000000 0.015625 0.000000000
[8,] 0.000000000 0.000000000 0.000000000 0.00000000 0.00000000 0.00000000 0.000000 0.009708738
>
> # inverse matrix of NRM
> ai = ainv(pedi)
>
> # print
> 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
>
> # LHS, RHS
>
> # LHS
> LHS = rbind(
+ cbind(t(X) %*% Rinv %*% X, t(X) %*% Rinv %*% Z, t(X) %*% Rinv %*% W),
+ cbind(t(Z) %*% Rinv %*% X, t(Z) %*% Rinv %*% Z, t(Z) %*% Rinv %*% W),
+ cbind(t(W) %*% Rinv %*% X, t(W) %*% Rinv %*% Z, t(W) %*% Rinv %*% W + ai * alpha))
>
> # print
> LHS
snp1 snp2 snp3
0.076267880 -0.0292324941 0.0165285648 0.0259434918 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000
snp1 -0.029232494 0.0294468531 -0.0028682259 -0.0205676752 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000
snp2 0.016528565 -0.0028682259 0.0222329297 -0.0049215685 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000
snp3 0.025943492 -0.0205676752 -0.0049215685 0.0385101610 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000
0.000000000 0.0000000000 0.0000000000 0.0000000000 10.428194 0.000000 3.476065 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000
0.000000000 0.0000000000 0.0000000000 0.0000000000 0.000000 10.428194 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000
0.000000000 0.0000000000 0.0000000000 0.0000000000 3.476065 0.000000 10.428194 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000
0.000000000 0.0000000000 0.0000000000 0.0000000000 0.000000 0.000000 0.000000 10.428194 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000
0.000000000 0.0000000000 0.0000000000 0.0000000000 0.000000 0.000000 0.000000 0.000000 10.428194 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000
0.000000000 0.0000000000 0.0000000000 0.0000000000 0.000000 0.000000 0.000000 0.000000 0.000000 10.428194 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000
0.000000000 0.0000000000 0.0000000000 0.0000000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 10.428194 0.000000 0.00000 0.000000 0.000000 0.000000
0.000000000 0.0000000000 0.0000000000 0.0000000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 10.428194 0.00000 0.000000 0.000000 0.000000
0.000000000 0.0000000000 0.0000000000 0.0000000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 13.90426 0.000000 0.000000 0.000000
0.000000000 0.0000000000 0.0000000000 0.0000000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 10.428194 0.000000 0.000000
0.000000000 0.0000000000 0.0000000000 0.0000000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 10.428194 0.000000
0.000000000 0.0000000000 0.0000000000 0.0000000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 10.428194
0.001792115 0.0024321557 -0.0006400410 0.0005120328 0.000000 0.000000 0.000000 3.476065 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000
0.001385042 0.0004946577 -0.0004946577 -0.0009893154 0.000000 0.000000 0.000000 0.000000 0.000000 3.476065 3.476065 3.476065 6.95213 3.476065 3.476065 3.476065
0.003333333 0.0011904762 0.0021428571 0.0042857143 0.000000 3.476065 0.000000 -6.952130 3.476065 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000
0.013698630 -0.0088062622 -0.0048923679 0.0176125245 0.000000 -6.952130 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000
0.019230769 -0.0123626374 0.0123626374 0.0054945055 0.000000 0.000000 0.000000 0.000000 -6.952130 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000
0.011494253 0.0041050903 0.0073891626 -0.0082101806 0.000000 0.000000 0.000000 0.000000 0.000000 -6.952130 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000
0.015625000 -0.0100446429 -0.0055803571 0.0044642857 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -6.95213 0.000000 0.000000 0.000000
0.009708738 -0.0062413315 0.0062413315 0.0027739251 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -6.95213 0.000000 0.000000 0.000000
0.000000000 0.0000000000 0.0000000000 0.0000000000 -6.952130 0.000000 -6.952130 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000
0.000000000 0.0000000000 0.0000000000 0.0000000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -6.952130 0.00000 0.000000 0.000000 0.000000
0.000000000 0.0000000000 0.0000000000 0.0000000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 -6.952130 0.000000
0.000000000 0.0000000000 0.0000000000 0.0000000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 -6.952130 0.000000 0.000000
0.000000000 0.0000000000 0.0000000000 0.0000000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -6.952130 0.000000 0.00000 0.000000 0.000000 0.000000
0.000000000 0.0000000000 0.0000000000 0.0000000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 -6.952130
0.0017921147 0.0013850416 0.003333333 0.013698630 0.019230769 0.011494253 0.015625000 0.009708738 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
snp1 0.0024321557 0.0004946577 0.001190476 -0.008806262 -0.012362637 0.004105090 -0.010044643 -0.006241331 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
snp2 -0.0006400410 -0.0004946577 0.002142857 -0.004892368 0.012362637 0.007389163 -0.005580357 0.006241331 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
snp3 0.0005120328 -0.0009893154 0.004285714 0.017612524 0.005494505 -0.008210181 0.004464286 0.002773925 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
0.0000000000 0.0000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 -6.95213 0.00000 0.00000 0.00000 0.00000 0.00000
0.0000000000 0.0000000000 3.476064811 -6.952129622 0.000000000 0.000000000 0.000000000 0.000000000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
0.0000000000 0.0000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 -6.95213 0.00000 0.00000 0.00000 0.00000 0.00000
3.4760648109 0.0000000000 -6.952129622 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
0.0000000000 0.0000000000 3.476064811 0.000000000 -6.952129622 0.000000000 0.000000000 0.000000000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
0.0000000000 3.4760648109 0.000000000 0.000000000 0.000000000 -6.952129622 0.000000000 0.000000000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
0.0000000000 3.4760648109 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.00000 0.00000 0.00000 0.00000 -6.95213 0.00000
0.0000000000 3.4760648109 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.00000 -6.95213 0.00000 0.00000 0.00000 0.00000
0.0000000000 6.9521296217 0.000000000 0.000000000 0.000000000 0.000000000 -6.952129622 -6.952129622 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
0.0000000000 3.4760648109 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.00000 0.00000 0.00000 -6.95213 0.00000 0.00000
0.0000000000 3.4760648109 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.00000 0.00000 -6.95213 0.00000 0.00000 0.00000
0.0000000000 3.4760648109 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.00000 0.00000 0.00000 0.00000 0.00000 -6.95213
10.4299865473 0.0000000000 -6.952129622 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
0.0000000000 34.7620331503 0.000000000 0.000000000 0.000000000 -6.952129622 -6.952129622 -6.952129622 0.00000 -6.95213 -6.95213 -6.95213 -6.95213 -6.95213
-6.9521296217 0.0000000000 20.859722199 -6.952129622 -6.952129622 0.000000000 0.000000000 0.000000000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
0.0000000000 0.0000000000 -6.952129622 13.917957874 0.000000000 0.000000000 0.000000000 0.000000000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
0.0000000000 0.0000000000 -6.952129622 0.000000000 13.923490013 0.000000000 0.000000000 0.000000000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
0.0000000000 -6.9521296217 0.000000000 0.000000000 0.000000000 13.915753496 0.000000000 0.000000000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
0.0000000000 -6.9521296217 0.000000000 0.000000000 0.000000000 0.000000000 13.919884243 0.000000000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
0.0000000000 -6.9521296217 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 13.913967981 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
0.0000000000 0.0000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 13.90426 0.00000 0.00000 0.00000 0.00000 0.00000
0.0000000000 -6.9521296217 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.00000 13.90426 0.00000 0.00000 0.00000 0.00000
0.0000000000 -6.9521296217 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.00000 0.00000 13.90426 0.00000 0.00000 0.00000
0.0000000000 -6.9521296217 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.00000 0.00000 0.00000 13.90426 0.00000 0.00000
0.0000000000 -6.9521296217 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.00000 0.00000 0.00000 0.00000 13.90426 0.00000
0.0000000000 -6.9521296217 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.00000 0.00000 0.00000 0.00000 0.00000 13.90426
>
> # RHS
> RHS = rbind(t(X) %*% Rinv %*% y,
+ t(Z) %*% Rinv %*% y,
+ t(W) %*% Rinv %*% y)
>
> # print
> RHS
[,1]
0.69592505
snp1 -0.26572369
snp2 0.04235790
snp3 0.34506266
0.00000000
0.00000000
0.00000000
0.00000000
0.00000000
0.00000000
0.00000000
0.00000000
0.00000000
0.00000000
0.00000000
0.00000000
0.01612903
0.01855956
0.04233333
0.21095890
0.11346154
0.08850575
0.15937500
0.04660194
0.00000000
0.00000000
0.00000000
0.00000000
0.00000000
0.00000000
>
> # Solutions
>
> # generalized inverse of LHS
> gi_LHS = ginv(LHS)
>
> # print
> gi_LHS
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,] 2.920047e+01 1.763942e+01 -2.231416e+01 -1.309408e+01 7.413843e-15 -8.839999e-03 -3.048551e-15 -5.289824e-03 3.231430e-04 -2.516429e-02 1.406606e-14 1.237420e-14
[2,] 1.763942e+01 6.933510e+01 1.431918e+00 2.532832e+01 -2.176580e-14 -5.022158e-03 -1.027103e-14 -7.728479e-03 2.593618e-02 -2.083596e-02 2.595898e-14 -1.328113e-14
[3,] -2.231416e+01 1.431918e+00 6.721137e+01 2.439554e+01 -1.212401e-14 1.563587e-02 -2.195249e-14 -2.340188e-02 -3.715804e-02 -3.300104e-03 -1.688944e-14 -3.007048e-14
[4,] -1.309408e+01 2.532832e+01 2.439554e+01 5.148543e+01 -3.281515e-14 -2.761546e-02 6.079561e-15 -3.312695e-02 -1.354714e-03 2.073423e-02 5.850616e-15 -2.586359e-14
[5,] 3.854270e-15 6.098944e-15 -1.627167e-14 3.412186e-15 1.438408e-01 3.281081e-17 1.665335e-16 2.460005e-17 1.294570e-18 4.900896e-17 -4.836998e-17 -7.003944e-17
[6,] -8.839999e-03 -5.022158e-03 1.563587e-02 -2.761546e-02 7.725859e-17 1.438160e-01 4.443646e-17 -1.373273e-06 4.761023e-06 -1.580866e-05 -7.508505e-17 -9.719437e-17
[7,] 2.268174e-15 -1.182284e-14 -7.722540e-15 -1.550289e-14 5.551115e-17 3.246030e-17 1.438408e-01 1.672703e-17 -5.036466e-19 -8.089353e-17 3.211128e-17 4.198234e-19
[8,] -5.289824e-03 -7.728479e-03 -2.340188e-02 -3.312695e-02 5.607399e-17 -1.373273e-06 2.890609e-17 1.438310e-01 -1.534947e-05 -4.746415e-07 -3.436971e-17 7.145836e-18
[9,] 3.231430e-04 2.593618e-02 -3.715804e-02 -1.354714e-03 -2.179995e-17 4.761023e-06 1.791440e-17 -1.534947e-05 1.437976e-01 1.099299e-05 -3.017503e-17 -1.800968e-17
[10,] -2.516429e-02 -2.083596e-02 -3.300104e-03 2.073423e-02 -4.667688e-17 -1.580866e-05 -1.394586e-18 -4.746415e-07 1.099299e-05 1.438223e-01 -4.292573e-17 1.362749e-17
[11,] 3.697444e-15 -1.750360e-16 -2.137471e-16 -1.108517e-14 -1.836047e-17 -5.174046e-17 -6.444924e-17 -1.799527e-17 -5.931322e-17 -1.385134e-16 1.438408e-01 -4.760164e-17
[12,] -1.280829e-14 -4.073019e-15 3.107244e-14 1.216392e-14 6.802225e-18 -6.472936e-17 7.657708e-17 -3.395347e-17 -4.383075e-17 8.026958e-18 1.631712e-16 1.438408e-01
[13,] -2.456054e-02 3.574275e-02 2.634502e-02 2.545956e-02 1.069449e-17 2.380274e-05 -8.101235e-18 1.889982e-05 3.217335e-05 1.078680e-05 -1.139568e-17 -4.671328e-17
[14,] -8.376406e-15 -1.839873e-14 1.950394e-14 -9.523788e-15 8.942065e-18 -2.578767e-17 8.775515e-17 -1.415649e-17 -4.335463e-17 4.526042e-18 8.900807e-17 -1.269771e-16
[15,] 4.888523e-15 2.162299e-14 9.652228e-15 1.121159e-14 -1.000579e-18 -8.186696e-17 3.638188e-17 -8.585356e-18 1.118782e-17 -2.005501e-17 5.235704e-17 -6.995308e-17
[16,] -8.825948e-16 -9.068977e-15 5.327753e-15 -9.894378e-15 2.158874e-17 5.648236e-18 -2.205240e-17 -1.871506e-17 -8.108127e-17 2.608850e-17 -2.733303e-16 1.974181e-16
[17,] -2.007289e-02 -3.825575e-02 -1.375616e-02 -4.014817e-02 3.817602e-18 6.134295e-06 3.610079e-17 -5.497659e-06 -2.782034e-05 1.146763e-05 2.945302e-17 5.202542e-17
[18,] -6.023642e-02 1.016341e-02 3.563529e-02 5.605151e-02 -2.221832e-17 7.293851e-06 6.619617e-19 1.364703e-05 3.841225e-05 1.517008e-06 -4.453104e-18 -1.423568e-16
[19,] -1.797118e-02 -3.072059e-02 -4.198090e-02 -6.976451e-02 5.844438e-17 1.007238e-06 3.244858e-17 7.190288e-02 -3.693437e-05 5.021854e-06 -1.637519e-17 3.630498e-17
[20,] -2.224559e-02 -2.289353e-02 2.463362e-03 -7.630544e-02 1.122973e-16 7.188370e-02 1.802412e-17 3.594938e-02 -1.132565e-05 -2.120206e-05 -4.135747e-17 -1.086394e-17
[21,] -8.500876e-03 2.354398e-02 -7.672751e-02 -3.691433e-02 1.154287e-17 7.645153e-06 4.597242e-17 3.592842e-02 7.183719e-02 1.900041e-05 -5.557656e-17 -4.033050e-17
[22,] -6.786464e-02 -2.617224e-02 1.286749e-02 5.912710e-02 -6.032017e-17 -2.006606e-05 2.865929e-17 6.111555e-06 3.569561e-05 7.189344e-02 -4.421854e-17 -5.923113e-17
[23,] -6.714570e-02 4.561056e-02 7.615665e-02 6.695275e-02 -1.666986e-17 3.695460e-05 -2.783898e-17 1.785869e-05 3.914425e-05 1.138376e-05 -1.404479e-17 -1.295412e-16
[24,] -4.221180e-02 3.603835e-02 1.216868e-02 4.001788e-02 2.115577e-17 1.794473e-05 -1.166689e-17 3.358797e-05 6.361469e-05 1.170686e-05 7.854100e-18 -9.209185e-17
[25,] 3.727857e-15 4.486850e-15 6.169878e-15 -7.616713e-15 7.192041e-02 7.935595e-17 7.192041e-02 2.590943e-17 -2.459203e-17 -1.294269e-17 3.537566e-17 -2.335057e-17
[26,] -3.011821e-02 5.081706e-03 1.781764e-02 2.802575e-02 3.327458e-19 3.646926e-06 3.540119e-17 6.823517e-06 1.920612e-05 7.585038e-07 1.180998e-16 7.192041e-02
[27,] -3.011821e-02 5.081706e-03 1.781764e-02 2.802575e-02 6.046083e-18 3.646926e-06 3.012527e-17 6.823517e-06 1.920612e-05 7.585038e-07 -2.587482e-17 -1.208657e-16
[28,] -3.011821e-02 5.081706e-03 1.781764e-02 2.802575e-02 -5.723745e-18 3.646926e-06 6.832831e-17 6.823517e-06 1.920612e-05 7.585038e-07 6.258824e-17 -1.330166e-16
[29,] -3.011821e-02 5.081706e-03 1.781764e-02 2.802575e-02 -5.791035e-17 3.646926e-06 -4.420883e-17 6.823517e-06 1.920612e-05 7.585038e-07 7.192041e-02 -1.312770e-16
[30,] -3.011821e-02 5.081706e-03 1.781764e-02 2.802575e-02 -6.157599e-18 3.646926e-06 -1.506816e-17 6.823517e-06 1.920612e-05 7.585038e-07 -1.837357e-16 1.964591e-17
[,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
[1,] -2.456054e-02 1.127386e-14 2.407494e-14 1.531248e-14 -2.007289e-02 -6.023642e-02 -1.797118e-02 -2.224559e-02 -8.500876e-03 -6.786464e-02 -6.714570e-02 -4.221180e-02
[2,] 3.574275e-02 3.171247e-14 7.965550e-15 -3.177764e-15 -3.825575e-02 1.016341e-02 -3.072059e-02 -2.289353e-02 2.354398e-02 -2.617224e-02 4.561056e-02 3.603835e-02
[3,] 2.634502e-02 -1.080795e-14 -4.995003e-14 -2.818954e-14 -1.375616e-02 3.563529e-02 -4.198090e-02 2.463362e-03 -7.672751e-02 1.286749e-02 7.615665e-02 1.216868e-02
[4,] 2.545956e-02 1.475052e-14 -2.039614e-14 -2.245477e-14 -4.014817e-02 5.605151e-02 -6.976451e-02 -7.630544e-02 -3.691433e-02 5.912710e-02 6.695275e-02 4.001788e-02
[5,] 7.721265e-17 1.005990e-17 3.408727e-17 6.035387e-18 3.991369e-17 3.686344e-17 5.348801e-17 3.946867e-17 5.312898e-17 6.283395e-17 3.211329e-17 7.978079e-17
[6,] 2.380274e-05 -3.892998e-17 -7.994268e-17 -5.936226e-17 6.134295e-06 7.293851e-06 1.007238e-06 7.188370e-02 7.645153e-06 -2.006606e-05 3.695460e-05 1.794473e-05
[7,] -2.265078e-17 2.104700e-17 -6.244156e-17 2.113832e-17 4.339012e-17 -4.803139e-17 4.450964e-17 4.576585e-17 2.170105e-17 -1.199753e-16 -2.480241e-17 -4.045046e-17
[8,] 1.889982e-05 -1.966336e-17 -1.981620e-18 1.267667e-17 -5.497659e-06 1.364703e-05 7.190288e-02 3.594938e-02 3.592842e-02 6.111555e-06 1.785869e-05 3.358797e-05
[9,] 3.217335e-05 -4.649104e-17 2.781690e-17 -2.191223e-17 -2.782034e-05 3.841225e-05 -3.693437e-05 -1.132565e-05 7.183719e-02 3.569561e-05 3.914425e-05 6.361469e-05
[10,] 1.078680e-05 2.953518e-17 3.697114e-18 2.689264e-18 1.146763e-05 1.517008e-06 5.021854e-06 -2.120206e-05 1.900041e-05 7.189344e-02 1.138376e-05 1.170686e-05
[11,] -4.621383e-17 -9.202127e-18 -1.000710e-16 -7.535890e-17 1.919274e-16 -1.258782e-16 1.246395e-16 8.698977e-17 -1.313935e-17 -1.451779e-16 -1.011369e-16 -8.747060e-17
[12,] 2.444268e-17 3.068831e-17 -1.176433e-16 -5.303253e-17 1.592858e-17 5.622370e-18 -2.659208e-17 -6.020184e-17 -5.736620e-17 2.424513e-17 2.499326e-17 1.070261e-17
[13,] 1.437820e-01 7.157615e-18 -5.193655e-17 -2.778831e-20 1.327353e-05 -4.014701e-05 3.498649e-05 5.319735e-05 6.575327e-05 -3.893299e-06 7.184602e-02 7.183707e-02
[14,] 9.835869e-18 1.438408e-01 2.676875e-17 -1.140285e-17 1.819425e-17 7.153431e-18 -5.566045e-18 -3.965451e-18 -5.302593e-17 1.067757e-17 6.415434e-18 6.785654e-18
[15,] -1.824879e-17 -1.073689e-17 1.438408e-01 1.461301e-17 1.314185e-18 9.103037e-18 -1.665990e-17 -8.136608e-17 1.034184e-17 -1.102141e-17 -3.141744e-18 -5.766190e-18
[16,] -6.346174e-18 -1.180660e-16 1.740054e-16 1.438408e-01 3.596347e-17 -4.121642e-17 7.420429e-18 1.633791e-17 -6.164964e-17 -9.038260e-18 -3.118813e-17 -2.160013e-17
[17,] 1.327353e-05 -3.644435e-18 2.772130e-18 2.535565e-17 1.438185e-01 2.476521e-05 7.190100e-02 3.595970e-02 3.590877e-02 2.958406e-05 2.128556e-05 3.002672e-05
[18,] -4.014701e-05 -1.672890e-20 -3.447643e-17 -6.134678e-17 2.476521e-05 1.437953e-01 3.285316e-05 2.736736e-05 7.404495e-05 7.189994e-02 7.186817e-02 7.184687e-02
[19,] 3.498649e-05 -3.113023e-17 5.688422e-18 2.234742e-17 7.190100e-02 3.285316e-05 1.438048e-01 7.190392e-02 7.184701e-02 2.395936e-05 3.743082e-05 6.539532e-05
[20,] 5.319735e-05 -5.675340e-17 -6.318892e-17 -2.293598e-17 3.595970e-02 2.736736e-05 7.190392e-02 1.437775e-01 3.593497e-02 -1.811941e-05 7.414731e-05 5.961475e-05
[21,] 6.575327e-05 -4.322887e-17 3.058577e-17 -5.851582e-18 3.590877e-02 7.404495e-05 7.184701e-02 3.593497e-02 1.436793e-01 6.552310e-05 7.743179e-05 1.281197e-04
[22,] -3.893299e-06 1.051766e-17 -2.497816e-17 -4.530919e-17 2.958406e-05 7.189994e-02 2.395936e-05 -1.811941e-05 6.552310e-05 1.437901e-01 3.595116e-02 3.594099e-02
[23,] 7.184602e-02 1.927940e-18 -8.488084e-17 -5.620609e-17 2.128556e-05 7.186817e-02 3.743082e-05 7.414731e-05 7.743179e-05 3.595116e-02 1.437335e-01 7.182676e-02
[24,] 7.183707e-02 1.346413e-17 -2.544736e-17 -1.621100e-17 3.002672e-05 7.184687e-02 6.539532e-05 5.961475e-05 1.281197e-04 3.594099e-02 7.182676e-02 1.436942e-01
[25,] 5.213887e-17 -8.519016e-18 -1.384489e-17 2.878586e-17 3.561158e-17 -5.831951e-18 4.669348e-17 8.592713e-17 1.235270e-17 -4.103907e-17 4.746855e-17 2.044665e-17
[26,] -2.007351e-05 3.359176e-17 -1.423844e-16 -2.746152e-17 1.238261e-05 7.189766e-02 1.642658e-05 1.368368e-05 3.702248e-05 3.594997e-02 3.593408e-02 3.592343e-02
[27,] -2.007351e-05 -1.150362e-17 7.192041e-02 -8.296526e-17 1.238261e-05 7.189766e-02 1.642658e-05 1.368368e-05 3.702248e-05 3.594997e-02 3.593408e-02 3.592343e-02
[28,] -2.007351e-05 7.192041e-02 -8.687369e-17 -1.184944e-17 1.238261e-05 7.189766e-02 1.642658e-05 1.368368e-05 3.702248e-05 3.594997e-02 3.593408e-02 3.592343e-02
[29,] -2.007351e-05 5.840608e-18 -8.513414e-17 -7.776410e-17 1.238261e-05 7.189766e-02 1.642658e-05 1.368368e-05 3.702248e-05 3.594997e-02 3.593408e-02 3.592343e-02
[30,] -2.007351e-05 -1.155881e-16 1.022180e-16 7.192041e-02 1.238261e-05 7.189766e-02 1.642658e-05 1.368368e-05 3.702248e-05 3.594997e-02 3.593408e-02 3.592343e-02
[,25] [,26] [,27] [,28] [,29] [,30]
[1,] 2.332834e-15 -3.011821e-02 -3.011821e-02 -3.011821e-02 -3.011821e-02 -3.011821e-02
[2,] -1.759614e-14 5.081706e-03 5.081706e-03 5.081706e-03 5.081706e-03 5.081706e-03
[3,] -1.645643e-14 1.781764e-02 1.781764e-02 1.781764e-02 1.781764e-02 1.781764e-02
[4,] -1.602331e-14 2.802575e-02 2.802575e-02 2.802575e-02 2.802575e-02 2.802575e-02
[5,] 7.192041e-02 -8.224665e-17 -5.685184e-17 4.062549e-17 1.599423e-16 1.556634e-17
[6,] 2.791043e-18 3.646926e-06 3.646926e-06 3.646926e-06 3.646926e-06 3.646926e-06
[7,] 7.192041e-02 4.094171e-17 6.546335e-17 6.389980e-18 -1.687216e-16 -5.119257e-17
[8,] 2.163392e-17 6.823517e-06 6.823517e-06 6.823517e-06 6.823517e-06 6.823517e-06
[9,] 3.497996e-17 1.920612e-05 1.920612e-05 1.920612e-05 1.920612e-05 1.920612e-05
[10,] -5.333961e-18 7.585038e-07 7.585038e-07 7.585038e-07 7.585038e-07 7.585038e-07
[11,] 1.417933e-16 -5.881716e-17 -1.334137e-16 -2.237553e-17 7.192041e-02 -1.420786e-16
[12,] 3.270198e-18 7.192041e-02 2.645242e-16 -6.852686e-17 -3.530294e-16 3.207997e-17
[13,] -5.204749e-20 -2.007351e-05 -2.007351e-05 -2.007351e-05 -2.007351e-05 -2.007351e-05
[14,] -8.856996e-17 -6.257526e-17 8.660753e-17 7.192041e-02 1.549181e-17 1.028841e-17
[15,] 3.298604e-17 -6.074041e-17 7.192041e-02 -3.991132e-17 3.987806e-17 1.385799e-17
[16,] -3.839556e-17 -2.951139e-17 -3.504387e-16 2.601217e-17 3.191725e-16 7.192041e-02
[17,] -3.192034e-18 1.238261e-05 1.238261e-05 1.238261e-05 1.238261e-05 1.238261e-05
[18,] -5.365654e-18 7.189766e-02 7.189766e-02 7.189766e-02 7.189766e-02 7.189766e-02
[19,] 1.009985e-17 1.642658e-05 1.642658e-05 1.642658e-05 1.642658e-05 1.642658e-05
[20,] 3.183542e-18 1.368368e-05 1.368368e-05 1.368368e-05 1.368368e-05 1.368368e-05
[21,] 3.382015e-17 3.702248e-05 3.702248e-05 3.702248e-05 3.702248e-05 3.702248e-05
[22,] 7.859526e-18 3.594997e-02 3.594997e-02 3.594997e-02 3.594997e-02 3.594997e-02
[23,] -3.735889e-17 3.593408e-02 3.593408e-02 3.593408e-02 3.593408e-02 3.593408e-02
[24,] 6.408456e-18 3.592343e-02 3.592343e-02 3.592343e-02 3.592343e-02 3.592343e-02
[25,] 1.438408e-01 -1.738295e-17 1.636175e-18 1.479296e-17 2.889214e-17 -2.276950e-17
[26,] -8.319072e-18 1.438294e-01 3.594883e-02 3.594883e-02 3.594883e-02 3.594883e-02
[27,] 3.244626e-17 3.594883e-02 1.438294e-01 3.594883e-02 3.594883e-02 3.594883e-02
[28,] -7.220066e-17 3.594883e-02 3.594883e-02 1.438294e-01 3.594883e-02 3.594883e-02
[29,] 9.163404e-17 3.594883e-02 3.594883e-02 3.594883e-02 1.438294e-01 3.594883e-02
[30,] -2.102841e-17 3.594883e-02 3.594883e-02 3.594883e-02 3.594883e-02 1.438294e-01
>
> # solution
> sol = gi_LHS %*% RHS
>
> # print
> sol
[,1]
[1,] 1.014413e+01
[2,] 2.654992e+00
[3,] -4.640234e+00
[4,] 2.951136e+00
[5,] 1.582177e-15
[6,] 1.486369e-03
[7,] -9.591387e-16
[8,] 6.595328e-04
[9,] -5.489412e-04
[10,] 1.402426e-03
[11,] -1.224625e-15
[12,] -2.331180e-15
[13,] -1.871431e-03
[14,] -3.404909e-15
[15,] 1.915600e-15
[16,] -1.403082e-15
[17,] -1.208474e-03
[18,] 8.051891e-05
[19,] 3.850622e-04
[20,] 2.422084e-03
[21,] -6.308807e-04
[22,] 2.143898e-03
[23,] -1.722584e-03
[24,] -1.939760e-03
[25,] -9.379962e-16
[26,] 4.025946e-05
[27,] 4.025946e-05
[28,] 4.025946e-05
[29,] 4.025946e-05
[30,] 4.025946e-05
>
> # mean effect of solutions
> mean_eff = sol[1]
> mean_eff
[1] 10.14413
>
> # snp effect of solutions
> snp_eff = sol[2:4]
> snp_eff
[1] 2.654992 -4.640234 2.951136
>
> # polygenic effect of solutions
> polge_eff = sol[5:30]
> polge_eff
[1] 1.582177e-15 1.486369e-03 -9.591387e-16 6.595328e-04 -5.489412e-04 1.402426e-03 -1.224625e-15 -2.331180e-15 -1.871431e-03 -3.404909e-15 1.915600e-15 -1.403082e-15
[13] -1.208474e-03 8.051891e-05 3.850622e-04 2.422084e-03 -6.308807e-04 2.143898e-03 -1.722584e-03 -1.939760e-03 -9.379962e-16 4.025946e-05 4.025946e-05 4.025946e-05
[25] 4.025946e-05 4.025946e-05
>
> # calculate DGV
> # genotype of animals
>
> # design matrix of animals
> P2 = sweep(M_all, 2, snp_mean * 2)
> P2
snp1 snp2 snp3
[1,] 1.3571429 -0.3571429 0.2857143
[2,] 0.3571429 -0.3571429 -0.7142857
[3,] 0.3571429 0.6428571 1.2857143
[4,] -0.6428571 -0.3571429 1.2857143
[5,] -0.6428571 0.6428571 0.2857143
[6,] 0.3571429 0.6428571 -0.7142857
[7,] -0.6428571 -0.3571429 0.2857143
[8,] -0.6428571 0.6428571 0.2857143
[9,] 1.3571429 -0.3571429 -0.7142857
[10,] -0.6428571 -0.3571429 -0.7142857
[11,] -0.6428571 0.6428571 0.2857143
[12,] 0.3571429 -0.3571429 -0.7142857
[13,] -0.6428571 -0.3571429 -0.7142857
[14,] 0.3571429 -0.3571429 0.2857143
>
> # DGV
> dgv = P2 %*% snp_eff
> dgv
[,1]
[1,] 6.1036121
[2,] 0.4974832
[3,] 1.7595226
[4,] 3.7447638
[5,] -3.8466063
[6,] -4.1427504
[7,] 0.7936273
[8,] -3.8466063
[9,] 3.1524756
[10,] -2.1575092
[11,] -3.8466063
[12,] 0.4974832
[13,] -2.1575092
[14,] 3.4486197
>
> # GEBV : animal 13 ~ 26
> gebv = dgv + polge_eff[13:26]
> gebv
[,1]
[1,] 6.1024036
[2,] 0.4975637
[3,] 1.7599076
[4,] 3.7471859
[5,] -3.8472371
[6,] -4.1406065
[7,] 0.7919048
[8,] -3.8485460
[9,] 3.1524756
[10,] -2.1574689
[11,] -3.8465660
[12,] 0.4975235
[13,] -2.1574689
[14,] 3.4486600
>
결과 책과 약간 다른다. Mean effect가 책과 결과가 다른데 책의 오타인 것으로 보인다. 그리고 DGV도 좀 다른데 round-off error로 보인다.
SNP effect를 고정효과로 처리하여 모형에 넣고, SNP effect를 추정한 다음 개체의 SNP effect의 design matrix 이용하여 개체의 DGV를 추정하였다.
'Animal Breeding > R for Genetic Evaluation' 카테고리의 다른 글
| Mixed Linear Model(SNP-BLUP Model) for Computing SNP effect, weighted analysis (0) | 2020.12.19 |
|---|---|
| Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects, unweighted analysis (0) | 2020.12.19 |
| Fixed Effect Model for SNP Effects, unweighted analysis (0) | 2020.12.16 |
| 임의 회귀 모형의 해를 이용한 개체 육종가 그래프 및 305일령 육종가 계산 (0) | 2020.11.26 |
| 다형질 모형의 차원 줄이기(Reduce the Dimension of Multivariate Models) (0) | 2020.09.29 |