# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.
# Raphael Mrode
# Example 11.1 p180
SNP effect를 고정효과로 모형에 넣어 SNP 추정하고 개체의 genotype을 이용하여 개체의 DGV(direct genomic breeding value)를 추정할 수 있다. SNP effect로 개체의 모든 유전적 효과를 설명할 수 없으므로 polygenic effect를 추가하여 분석한다. 그래서 분석모형은 y = mean + snp effect + polygenic effect + e 가 된다.
SNP effect의 design matrix를 만드는 방법은 다음과 같다. SNP gentype을 먼저 0, 1, 2로 코드화 한다. SNP 별로 평균을 구하는데 각 개체별로 두 개의 SNP가 동원된 것으므로 SNP 별 합을 개체의 2배수로 나누어 평균을 구한다. 코드화한 genotype에서 각 SNP의 평균의 2배를 빼 SNP effect의 design matrix를 구한다. 이 개체별 design matrix를 이용하여 나중에 SNP effect 추정 후 DGV를 추정한다.
예제에서는 단지 3개의 SNP만을 이용하는 것을 보여준다.
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, unweighted 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
# inverse matrix of NRM
ai = ainv(pedi)
# print
ai
# LHS, RHS
# LHS
LHS = rbind(
cbind(t(X) %*% X, t(X) %*% Z, t(X) %*% W),
cbind(t(Z) %*% X, t(Z) %*% Z, t(Z) %*% W),
cbind(t(W) %*% X, t(W) %*% Z, t(W) %*% W + ai * alpha))
# print
LHS
# RHS
RHS = rbind(t(X) %*% y,
t(Z) %*% y,
t(W) %*% 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, unweighted 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/11_fixed effecct model for SNP effects_unweighted"
>
> # 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
>
> # 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]
[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
[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
[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
[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
[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
[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
[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
[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
[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
[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
[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
[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
[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
[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
[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
[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
[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
[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
[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
[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
[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
[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
[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
[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
[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
[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
[,26]
[1,] 0
[2,] 0
[3,] 0
[4,] 0
[5,] 0
[6,] 0
[7,] 0
[8,] 0
[9,] 0
[10,] 0
[11,] 0
[12,] -1
[13,] 0
[14,] -1
[15,] 0
[16,] 0
[17,] 0
[18,] 0
[19,] 0
[20,] 0
[21,] 0
[22,] 0
[23,] 0
[24,] 0
[25,] 0
[26,] 2
>
> # LHS, RHS
>
> # LHS
> LHS = rbind(
+ cbind(t(X) %*% X, t(X) %*% Z, t(X) %*% W),
+ cbind(t(Z) %*% X, t(Z) %*% Z, t(Z) %*% W),
+ cbind(t(W) %*% X, t(W) %*% Z, t(W) %*% W + ai * alpha))
>
> # print
> LHS
snp1 snp2 snp3
8.0000000 -0.1428571 1.1428571 2.2857143 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000
snp1 -0.1428571 3.8775510 -0.5204082 -1.0408163 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000
snp2 1.1428571 -0.5204082 2.1632653 0.3265306 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000
snp3 2.2857143 -1.0408163 0.3265306 4.6530612 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000
0.0000000 0.0000000 0.0000000 0.0000000 10.428194 0.000000 3.476065 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000
0.0000000 0.0000000 0.0000000 0.0000000 0.000000 10.428194 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000
0.0000000 0.0000000 0.0000000 0.0000000 3.476065 0.000000 10.428194 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000
0.0000000 0.0000000 0.0000000 0.0000000 0.000000 0.000000 0.000000 10.428194 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000
0.0000000 0.0000000 0.0000000 0.0000000 0.000000 0.000000 0.000000 0.000000 10.428194 0.000000 0.000000 0.000000 0.00000 0.000000
0.0000000 0.0000000 0.0000000 0.0000000 0.000000 0.000000 0.000000 0.000000 0.000000 10.428194 0.000000 0.000000 0.00000 0.000000
0.0000000 0.0000000 0.0000000 0.0000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 10.428194 0.000000 0.00000 0.000000
0.0000000 0.0000000 0.0000000 0.0000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 10.428194 0.00000 0.000000
0.0000000 0.0000000 0.0000000 0.0000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 13.90426 0.000000
0.0000000 0.0000000 0.0000000 0.0000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 10.428194
0.0000000 0.0000000 0.0000000 0.0000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000
0.0000000 0.0000000 0.0000000 0.0000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000
1.0000000 1.3571429 -0.3571429 0.2857143 0.000000 0.000000 0.000000 3.476065 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000
1.0000000 0.3571429 -0.3571429 -0.7142857 0.000000 0.000000 0.000000 0.000000 0.000000 3.476065 3.476065 3.476065 6.95213 3.476065
1.0000000 0.3571429 0.6428571 1.2857143 0.000000 3.476065 0.000000 -6.952130 3.476065 0.000000 0.000000 0.000000 0.00000 0.000000
1.0000000 -0.6428571 -0.3571429 1.2857143 0.000000 -6.952130 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000
1.0000000 -0.6428571 0.6428571 0.2857143 0.000000 0.000000 0.000000 0.000000 -6.952130 0.000000 0.000000 0.000000 0.00000 0.000000
1.0000000 0.3571429 0.6428571 -0.7142857 0.000000 0.000000 0.000000 0.000000 0.000000 -6.952130 0.000000 0.000000 0.00000 0.000000
1.0000000 -0.6428571 -0.3571429 0.2857143 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -6.95213 0.000000
1.0000000 -0.6428571 0.6428571 0.2857143 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -6.95213 0.000000
0.0000000 0.0000000 0.0000000 0.0000000 -6.952130 0.000000 -6.952130 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000
0.0000000 0.0000000 0.0000000 0.0000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -6.952130 0.00000 0.000000
0.0000000 0.0000000 0.0000000 0.0000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000
0.0000000 0.0000000 0.0000000 0.0000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 -6.952130
0.0000000 0.0000000 0.0000000 0.0000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -6.952130 0.000000 0.00000 0.000000
0.0000000 0.0000000 0.0000000 0.0000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000
0.000000 0.000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 0.00000 0.00000 0.00000 0.00000
snp1 0.000000 0.000000 1.3571429 0.3571429 0.3571429 -0.6428571 -0.6428571 0.3571429 -0.6428571 -0.6428571 0.00000 0.00000 0.00000 0.00000
snp2 0.000000 0.000000 -0.3571429 -0.3571429 0.6428571 -0.3571429 0.6428571 0.6428571 -0.3571429 0.6428571 0.00000 0.00000 0.00000 0.00000
snp3 0.000000 0.000000 0.2857143 -0.7142857 1.2857143 1.2857143 0.2857143 -0.7142857 0.2857143 0.2857143 0.00000 0.00000 0.00000 0.00000
0.000000 0.000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 -6.95213 0.00000 0.00000 0.00000
0.000000 0.000000 0.0000000 0.0000000 3.4760648 -6.9521296 0.0000000 0.0000000 0.0000000 0.0000000 0.00000 0.00000 0.00000 0.00000
0.000000 0.000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 -6.95213 0.00000 0.00000 0.00000
0.000000 0.000000 3.4760648 0.0000000 -6.9521296 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00000 0.00000 0.00000 0.00000
0.000000 0.000000 0.0000000 0.0000000 3.4760648 0.0000000 -6.9521296 0.0000000 0.0000000 0.0000000 0.00000 0.00000 0.00000 0.00000
0.000000 0.000000 0.0000000 3.4760648 0.0000000 0.0000000 0.0000000 -6.9521296 0.0000000 0.0000000 0.00000 0.00000 0.00000 0.00000
0.000000 0.000000 0.0000000 3.4760648 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00000 0.00000 0.00000 0.00000
0.000000 0.000000 0.0000000 3.4760648 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00000 -6.95213 0.00000 0.00000
0.000000 0.000000 0.0000000 6.9521296 0.0000000 0.0000000 0.0000000 0.0000000 -6.9521296 -6.9521296 0.00000 0.00000 0.00000 0.00000
0.000000 0.000000 0.0000000 3.4760648 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00000 0.00000 0.00000 -6.95213
10.428194 0.000000 0.0000000 3.4760648 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00000 0.00000 -6.95213 0.00000
0.000000 10.428194 0.0000000 3.4760648 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00000 0.00000 0.00000 0.00000
0.000000 0.000000 11.4281944 0.0000000 -6.9521296 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00000 0.00000 0.00000 0.00000
3.476065 3.476065 0.0000000 35.7606481 0.0000000 0.0000000 0.0000000 -6.9521296 -6.9521296 -6.9521296 0.00000 -6.95213 -6.95213 -6.95213
0.000000 0.000000 -6.9521296 0.0000000 21.8563889 -6.9521296 -6.9521296 0.0000000 0.0000000 0.0000000 0.00000 0.00000 0.00000 0.00000
0.000000 0.000000 0.0000000 0.0000000 -6.9521296 14.9042592 0.0000000 0.0000000 0.0000000 0.0000000 0.00000 0.00000 0.00000 0.00000
0.000000 0.000000 0.0000000 0.0000000 -6.9521296 0.0000000 14.9042592 0.0000000 0.0000000 0.0000000 0.00000 0.00000 0.00000 0.00000
0.000000 0.000000 0.0000000 -6.9521296 0.0000000 0.0000000 0.0000000 14.9042592 0.0000000 0.0000000 0.00000 0.00000 0.00000 0.00000
0.000000 0.000000 0.0000000 -6.9521296 0.0000000 0.0000000 0.0000000 0.0000000 14.9042592 0.0000000 0.00000 0.00000 0.00000 0.00000
0.000000 0.000000 0.0000000 -6.9521296 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 14.9042592 0.00000 0.00000 0.00000 0.00000
0.000000 0.000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 13.90426 0.00000 0.00000 0.00000
0.000000 0.000000 0.0000000 -6.9521296 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00000 13.90426 0.00000 0.00000
-6.952130 0.000000 0.0000000 -6.9521296 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00000 0.00000 13.90426 0.00000
0.000000 0.000000 0.0000000 -6.9521296 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00000 0.00000 0.00000 13.90426
0.000000 0.000000 0.0000000 -6.9521296 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00000 0.00000 0.00000 0.00000
0.000000 -6.952130 0.0000000 -6.9521296 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.00000 0.00000 0.00000 0.00000
0.00000 0.00000
snp1 0.00000 0.00000
snp2 0.00000 0.00000
snp3 0.00000 0.00000
0.00000 0.00000
0.00000 0.00000
0.00000 0.00000
0.00000 0.00000
0.00000 0.00000
0.00000 0.00000
-6.95213 0.00000
0.00000 0.00000
0.00000 0.00000
0.00000 0.00000
0.00000 0.00000
0.00000 -6.95213
0.00000 0.00000
-6.95213 -6.95213
0.00000 0.00000
0.00000 0.00000
0.00000 0.00000
0.00000 0.00000
0.00000 0.00000
0.00000 0.00000
0.00000 0.00000
0.00000 0.00000
0.00000 0.00000
0.00000 0.00000
13.90426 0.00000
0.00000 13.90426
>
> # RHS
> RHS = rbind(t(X) %*% y,
+ t(Z) %*% y,
+ t(W) %*% y)
>
> # print
> RHS
[,1]
79.10
snp1 0.95
snp2 2.85
snp3 29.60
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
9.00
13.40
12.70
15.40
5.90
7.70
10.20
4.80
0.00
0.00
0.00
0.00
0.00
0.00
>
> # Solutions
>
> # generalized inverse of LHS
> gi_LHS = ginv(LHS)
>
> # print
> gi_LHS
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 2.029219e-01 -3.255082e-02 -8.500123e-02 -9.538670e-02 -8.813991e-17 -7.749019e-03 -5.907393e-17 -7.339499e-03 -8.423214e-03 -1.124651e-02
[2,] -3.255082e-02 3.248575e-01 8.134068e-02 9.542976e-02 -1.098782e-16 1.072783e-02 1.818259e-16 -4.884210e-03 1.164245e-02 -4.992010e-03
[3,] -8.500123e-02 8.134068e-02 5.604372e-01 2.307770e-02 4.164755e-16 2.081808e-02 1.597493e-16 -1.878828e-02 -1.406575e-02 -1.939051e-02
[4,] -9.538670e-02 9.542976e-02 2.307770e-02 3.295430e-01 -5.230087e-17 -1.481457e-02 5.306137e-17 -2.760843e-02 5.080273e-03 1.647388e-02
[5,] 2.329623e-18 -3.804909e-16 6.298383e-19 -2.879471e-16 1.438408e-01 4.622486e-16 1.249001e-16 -5.883868e-16 7.059733e-16 -2.124216e-18
[6,] -7.749019e-03 1.072783e-02 2.081808e-02 -1.481457e-02 1.074399e-17 1.418405e-01 5.505648e-17 -1.363038e-04 3.553867e-04 -1.322611e-03
[7,] -9.458291e-17 -2.158422e-16 1.088446e-16 -2.082423e-16 1.249001e-16 4.318604e-16 1.438408e-01 -6.173507e-16 7.389172e-16 -6.170158e-17
[8,] -7.339499e-03 -4.884210e-03 -1.878828e-02 -2.760843e-02 -6.340412e-17 -1.363038e-04 1.479604e-17 1.424464e-01 -7.088821e-04 1.078415e-04
[9,] -8.423214e-03 1.164245e-02 -1.406575e-02 5.080273e-03 2.970513e-17 3.553867e-04 3.015109e-17 -7.088821e-04 1.407702e-01 1.005601e-03
[10,] -1.124651e-02 -4.992010e-03 -1.939051e-02 1.647388e-02 -6.298026e-17 -1.322611e-03 7.610779e-18 1.078415e-04 1.005601e-03 1.415624e-01
[11,] 2.602085e-18 -8.933826e-17 -2.775558e-17 -2.263814e-16 -1.173397e-16 2.873136e-18 -1.536369e-16 6.017322e-17 -1.454864e-17 5.895349e-17
[12,] -2.862294e-17 2.602085e-17 2.428613e-17 -5.030698e-17 6.117076e-17 1.610040e-17 1.851744e-16 4.044074e-17 -3.342731e-17 1.669671e-17
[13,] -1.922866e-02 2.429536e-02 5.003082e-03 3.922659e-03 -7.208471e-17 1.958920e-03 -5.429396e-17 1.885193e-03 1.888782e-03 7.119691e-04
[14,] -1.214306e-17 6.765422e-17 -1.474515e-16 -4.336809e-18 -1.062664e-17 1.512462e-17 -1.474799e-16 4.380177e-17 5.109303e-18 -2.287667e-17
[15,] 3.469447e-18 -1.170938e-16 1.214306e-17 -9.020562e-17 1.516989e-16 -5.025277e-17 1.248652e-16 1.647987e-17 1.526015e-17 4.526544e-18
[16,] -1.821460e-17 7.546047e-17 -3.816392e-17 1.040834e-17 -4.032172e-17 -3.680866e-17 -3.779099e-18 7.166576e-17 -1.586323e-17 3.686287e-18
[17,] -2.676957e-02 -5.539375e-02 4.722536e-03 -3.923171e-02 1.830766e-16 4.916905e-04 -2.540570e-17 6.855658e-04 -2.361692e-03 8.977596e-04
[18,] -6.308434e-02 1.860433e-02 2.170085e-02 5.617790e-02 -6.347784e-17 6.532274e-04 3.333540e-17 -4.389663e-04 2.891378e-03 8.778743e-04
[19,] -2.439404e-02 -3.502319e-02 -2.582116e-02 -6.102850e-02 8.483571e-17 4.138959e-05 2.492510e-17 7.017152e-02 -2.244169e-03 6.106421e-04
[20,] -2.382055e-02 -1.419855e-03 1.831654e-02 -5.273610e-02 7.677376e-17 6.894064e-02 7.392803e-17 3.488130e-02 -5.890045e-04 -1.678596e-03
[21,] -2.483184e-02 -4.792207e-05 -3.400920e-02 -2.289384e-02 1.050041e-17 5.537749e-04 6.646180e-17 3.402244e-02 6.619246e-02 1.813723e-03
[22,] -4.841193e-02 1.814151e-03 -1.823534e-02 5.279976e-02 -7.729679e-17 -1.657303e-03 -3.493507e-18 -5.772085e-05 2.954091e-03 6.894169e-02
[23,] -5.362241e-02 3.632630e-02 3.465475e-02 3.278580e-02 -7.449488e-17 2.983928e-03 1.998763e-17 1.035410e-03 2.862601e-03 5.004040e-04
[24,] -4.791926e-02 3.086876e-02 -2.947737e-03 3.123741e-02 -1.102180e-16 1.587140e-03 -3.077451e-17 2.296009e-03 3.806341e-03 1.801408e-03
[25,] 7.898705e-18 -5.436645e-16 -2.488422e-17 -3.772473e-16 7.192041e-02 6.303627e-16 7.192041e-02 -8.573219e-16 1.034442e-15 -2.944094e-17
[26,] -3.154217e-02 9.302166e-03 1.085042e-02 2.808895e-02 2.137914e-17 3.266137e-04 1.717343e-16 -2.194831e-04 1.445689e-03 4.389372e-04
[27,] -3.154217e-02 9.302166e-03 1.085042e-02 2.808895e-02 5.746220e-17 3.266137e-04 1.583557e-16 -2.194831e-04 1.445689e-03 4.389372e-04
[28,] -3.154217e-02 9.302166e-03 1.085042e-02 2.808895e-02 -4.739571e-17 3.266137e-04 -8.756498e-17 -2.194831e-04 1.445689e-03 4.389372e-04
[29,] -3.154217e-02 9.302166e-03 1.085042e-02 2.808895e-02 -1.313957e-16 3.266137e-04 -9.889781e-17 -2.194831e-04 1.445689e-03 4.389372e-04
[30,] -3.154217e-02 9.302166e-03 1.085042e-02 2.808895e-02 -6.701577e-17 3.266137e-04 -2.937504e-17 -2.194831e-04 1.445689e-03 4.389372e-04
[,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
[1,] -5.551115e-17 -7.806256e-18 -1.922866e-02 5.204170e-18 -2.775558e-17 -2.515349e-17 -2.676957e-02 -6.308434e-02 -2.439404e-02 -2.382055e-02
[2,] -9.540979e-18 6.938894e-18 2.429536e-02 -8.760354e-17 1.066855e-16 -9.714451e-17 -5.539375e-02 1.860433e-02 -3.502319e-02 -1.419855e-03
[3,] -1.439820e-16 -1.561251e-17 5.003082e-03 -1.422473e-16 7.459311e-17 -1.214306e-17 4.722536e-03 2.170085e-02 -2.582116e-02 1.831654e-02
[4,] -1.821460e-17 -5.204170e-18 3.922659e-03 -9.194034e-17 3.729655e-17 -4.597017e-17 -3.923171e-02 5.617790e-02 -6.102850e-02 -5.273610e-02
[5,] -5.247009e-17 -1.715734e-16 2.485646e-16 -1.381342e-16 -1.409332e-16 -2.060732e-16 -7.631449e-17 8.574258e-17 -3.209790e-16 3.571539e-17
[6,] -2.883978e-17 -1.149254e-17 1.958920e-03 7.155734e-18 6.619054e-17 -3.393553e-17 4.916905e-04 6.532274e-04 4.138959e-05 6.894064e-02
[7,] -1.441928e-16 -8.548200e-17 3.535717e-16 -2.256486e-16 -2.257681e-16 -3.538410e-17 -2.422365e-17 2.345273e-16 -3.973950e-16 5.911108e-17
[8,] 4.857226e-17 5.724587e-17 1.885193e-03 6.071532e-17 2.320193e-17 5.822166e-17 6.855658e-04 -4.389663e-04 7.017152e-02 3.488130e-02
[9,] -7.264155e-18 3.416592e-17 1.888782e-03 -1.091656e-17 4.254138e-17 1.465706e-17 -2.361692e-03 2.891378e-03 -2.244169e-03 -5.890045e-04
[10,] 1.620882e-17 1.859407e-17 7.119691e-04 1.493488e-17 -1.414884e-17 3.930233e-18 8.977596e-04 8.778743e-04 6.106421e-04 -1.678596e-03
[11,] 1.438408e-01 -1.816039e-17 -2.753874e-17 3.339343e-17 -2.152141e-17 -7.746625e-17 1.052760e-16 1.604619e-17 1.476683e-16 9.432559e-17
[12,] 9.698188e-17 1.438408e-01 2.883978e-17 -8.164042e-17 2.222614e-18 -2.423192e-17 3.024924e-17 -3.903128e-18 6.526897e-17 6.112190e-17
[13,] 9.671083e-17 -6.635317e-17 1.390523e-01 3.512815e-17 -4.380177e-17 4.185020e-17 3.589691e-06 -1.659950e-03 2.829584e-03 4.353172e-03
[14,] -5.177065e-17 -5.187907e-17 5.312591e-17 1.438408e-01 -8.407988e-17 -1.859407e-17 -1.561251e-17 1.084202e-17 4.228388e-17 2.556007e-17
[15,] -6.911789e-17 -1.414884e-17 5.204170e-17 1.778092e-17 1.438408e-01 -6.559423e-18 3.469447e-17 -2.645453e-17 4.401861e-17 -8.782038e-18
[16,] 4.591596e-17 1.263096e-17 2.688821e-17 -1.530893e-16 -5.144539e-17 1.438408e-01 -8.565197e-18 -2.038300e-17 5.009014e-17 2.425902e-17
[17,] 3.946496e-17 2.949030e-17 3.589691e-06 4.130810e-17 -1.702197e-17 3.144186e-17 1.407936e-01 3.330344e-03 7.142513e-02 3.645010e-02
[18,] 1.023487e-16 -2.168404e-17 -1.659950e-03 3.122502e-17 6.505213e-18 2.645453e-17 3.330344e-03 1.381869e-01 1.006723e-03 1.483202e-03
[19,] 5.052382e-17 4.727121e-17 2.829584e-03 6.028164e-17 -6.071532e-18 4.531965e-17 7.142513e-02 1.006723e-03 1.409698e-01 7.054701e-02
[20,] 8.104411e-18 5.204170e-17 4.353172e-03 3.932943e-17 6.114900e-17 6.423898e-18 3.645010e-02 1.483202e-03 7.054701e-02 1.386845e-01
[21,] 2.580401e-17 1.886512e-17 4.247966e-03 1.951564e-17 8.890458e-18 2.092510e-17 3.217003e-02 4.840428e-03 6.711867e-02 3.439000e-02
[22,] 6.071532e-17 -8.239937e-18 2.379786e-04 2.298509e-17 -1.257675e-17 6.505213e-18 3.011811e-03 7.041027e-02 1.419324e-03 -1.776292e-03
[23,] 1.227317e-16 -6.461845e-17 6.646977e-02 2.688821e-17 -2.081668e-17 3.209238e-17 1.827191e-03 6.816151e-02 2.466711e-03 5.709247e-03
[24,] 1.079865e-16 -5.854692e-17 6.613409e-02 3.426079e-17 -3.035766e-17 3.339343e-17 1.510333e-03 6.670550e-02 4.199179e-03 4.480300e-03
[25,] -1.427243e-16 -1.969319e-16 3.742461e-16 -2.425523e-16 -2.478455e-16 -1.864753e-16 -5.631687e-17 1.414053e-16 -5.271757e-16 6.804336e-17
[26,] 1.951564e-17 7.192041e-02 -8.299750e-04 1.461505e-16 -1.821460e-17 -5.854692e-17 1.665172e-03 6.909345e-02 5.033613e-04 7.416012e-04
[27,] 9.150666e-17 -3.035766e-17 -8.299750e-04 -6.158268e-17 7.192041e-02 6.114900e-17 1.665172e-03 6.909345e-02 5.033613e-04 7.416012e-04
[28,] 6.028164e-17 -1.387779e-17 -8.299750e-04 7.192041e-02 -1.647987e-17 1.621966e-16 1.665172e-03 6.909345e-02 5.033613e-04 7.416012e-04
[29,] 7.192041e-02 -6.895526e-17 -8.299750e-04 -9.974660e-18 5.247539e-17 -6.938894e-17 1.665172e-03 6.909345e-02 5.033613e-04 7.416012e-04
[30,] -1.734723e-18 9.107298e-18 -8.299750e-04 1.162265e-16 -6.808790e-17 7.192041e-02 1.665172e-03 6.909345e-02 5.033613e-04 7.416012e-04
[,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30]
[1,] -2.483184e-02 -4.841193e-02 -5.362241e-02 -4.791926e-02 -1.050221e-16 -3.154217e-02 -3.154217e-02 -3.154217e-02 -3.154217e-02 -3.154217e-02
[2,] -4.792207e-05 1.814151e-03 3.632630e-02 3.086876e-02 9.667407e-17 9.302166e-03 9.302166e-03 9.302166e-03 9.302166e-03 9.302166e-03
[3,] -3.400920e-02 -1.823534e-02 3.465475e-02 -2.947737e-03 4.748412e-16 1.085042e-02 1.085042e-02 1.085042e-02 1.085042e-02 1.085042e-02
[4,] -2.289384e-02 5.279976e-02 3.278580e-02 3.123741e-02 4.566415e-18 2.808895e-02 2.808895e-02 2.808895e-02 2.808895e-02 2.808895e-02
[5,] 2.611025e-16 7.438382e-17 2.170978e-16 2.341855e-16 7.192041e-02 2.074505e-17 -1.276276e-16 -1.703319e-16 -3.906911e-17 2.854173e-17
[6,] 5.537749e-04 -1.657303e-03 2.983928e-03 1.587140e-03 6.913168e-17 3.266137e-04 3.266137e-04 3.266137e-04 3.266137e-04 3.266137e-04
[7,] 2.785272e-16 1.014383e-16 3.715807e-16 3.646591e-16 7.192041e-02 -1.650363e-17 -1.828810e-18 1.472670e-16 8.431150e-17 -1.323675e-17
[8,] 3.402244e-02 -5.772085e-05 1.035410e-03 2.296009e-03 2.494461e-17 -2.194831e-04 -2.194831e-04 -2.194831e-04 -2.194831e-04 -2.194831e-04
[9,] 6.619246e-02 2.954091e-03 2.862601e-03 3.806341e-03 5.572096e-17 1.445689e-03 1.445689e-03 1.445689e-03 1.445689e-03 1.445689e-03
[10,] 1.813723e-03 6.894169e-02 5.004040e-04 1.801408e-03 -4.553002e-17 4.389372e-04 4.389372e-04 4.389372e-04 4.389372e-04 4.389372e-04
[11,] 7.686993e-17 3.079134e-17 -4.770490e-18 -4.033232e-17 2.049676e-16 1.214306e-17 1.127570e-17 1.778092e-17 7.192041e-02 -6.765422e-17
[12,] 1.832302e-17 -1.474515e-17 9.974660e-18 -1.734723e-17 -8.632881e-17 7.192041e-02 1.040834e-17 -4.380177e-17 2.558717e-17 -2.341877e-17
[13,] 4.247966e-03 2.379786e-04 6.646977e-02 6.613409e-02 -8.332182e-17 -8.299750e-04 -8.299750e-04 -8.299750e-04 -8.299750e-04 -8.299750e-04
[14,] 4.976488e-17 -1.561251e-17 3.512815e-17 1.994932e-17 7.420757e-17 -4.336809e-17 -7.979728e-17 7.192041e-02 -3.946496e-17 2.255141e-17
[15,] 3.491131e-17 -3.165870e-17 7.372575e-18 -1.778092e-17 -2.976979e-16 -1.647987e-17 7.192041e-02 1.301043e-17 -5.854692e-17 -4.336809e-17
[16,] 1.593777e-17 -2.688821e-17 2.168404e-18 -1.257675e-17 1.265510e-16 -2.168404e-18 -2.255141e-17 -1.543904e-16 3.035766e-17 7.192041e-02
[17,] 3.217003e-02 3.011811e-03 1.827191e-03 1.510333e-03 5.811441e-17 1.665172e-03 1.665172e-03 1.665172e-03 1.665172e-03 1.665172e-03
[18,] 4.840428e-03 7.041027e-02 6.816151e-02 6.670550e-02 -1.879806e-17 6.909345e-02 6.909345e-02 6.909345e-02 6.909345e-02 6.909345e-02
[19,] 6.711867e-02 1.419324e-03 2.466711e-03 4.199179e-03 9.221979e-17 5.033613e-04 5.033613e-04 5.033613e-04 5.033613e-04 5.033613e-04
[20,] 3.439000e-02 -1.776292e-03 5.709247e-03 4.480300e-03 1.174630e-16 7.416012e-04 7.416012e-04 7.416012e-04 7.416012e-04 7.416012e-04
[21,] 1.328480e-01 5.140798e-03 5.527257e-03 7.809102e-03 6.274527e-17 2.420214e-03 2.420214e-03 2.420214e-03 2.420214e-03 2.420214e-03
[22,] 5.140798e-03 1.386177e-01 3.483136e-02 3.605486e-02 -7.342940e-17 3.520513e-02 3.520513e-02 3.520513e-02 3.520513e-02 3.520513e-02
[23,] 5.527257e-03 3.483136e-02 1.352606e-01 6.584049e-02 -3.876275e-17 3.408076e-02 3.408076e-02 3.408076e-02 3.408076e-02 3.408076e-02
[24,] 7.809102e-03 3.605486e-02 6.584049e-02 1.331332e-01 -8.314716e-17 3.335275e-02 3.335275e-02 3.335275e-02 3.335275e-02 3.335275e-02
[25,] 3.845722e-16 8.174885e-17 3.219778e-16 3.378198e-16 1.438408e-01 -5.498276e-17 -1.277190e-16 -4.527046e-17 -1.257900e-17 -3.813977e-17
[26,] 2.420214e-03 3.520513e-02 3.408076e-02 3.335275e-02 -6.326058e-17 1.424273e-01 3.454673e-02 3.454673e-02 3.454673e-02 3.454673e-02
[27,] 2.420214e-03 3.520513e-02 3.408076e-02 3.335275e-02 -2.253800e-16 3.454673e-02 1.424273e-01 3.454673e-02 3.454673e-02 3.454673e-02
[28,] 2.420214e-03 3.520513e-02 3.408076e-02 3.335275e-02 4.798244e-17 3.454673e-02 3.454673e-02 1.424273e-01 3.454673e-02 3.454673e-02
[29,] 2.420214e-03 3.520513e-02 3.408076e-02 3.335275e-02 1.484779e-16 3.454673e-02 3.454673e-02 3.454673e-02 1.424273e-01 3.454673e-02
[30,] 2.420214e-03 3.520513e-02 3.408076e-02 3.335275e-02 5.835403e-17 3.454673e-02 3.454673e-02 3.454673e-02 3.454673e-02 1.424273e-01
>
> # solution
> sol = gi_LHS %*% RHS
>
> # print
> sol
[,1]
[1,] 9.895358e+00
[2,] 6.068652e-01
[3,] -4.080274e+00
[4,] 1.934155e+00
[5,] -6.311179e-15
[6,] 1.220146e-01
[7,] -6.787405e-15
[8,] 1.945577e-01
[9,] -1.042586e-01
[10,] 9.507380e-02
[11,] -1.720163e-15
[12,] -1.656048e-15
[13,] -2.644430e-01
[14,] 1.180100e-16
[15,] -2.138676e-15
[16,] -6.615097e-16
[17,] -2.988163e-01
[18,] 2.558718e-01
[19,] 1.424285e-01
[20,] 2.542361e-01
[21,] -8.517367e-02
[22,] 2.705466e-01
[23,] -9.223368e-02
[24,] -1.807806e-01
[25,] -7.584268e-15
[26,] 1.279359e-01
[27,] 1.279359e-01
[28,] 1.279359e-01
[29,] 1.279359e-01
[30,] 1.279359e-01
>
> # mean effect of solutions
> mean_eff = sol[1]
> mean_eff
[1] 9.895358
>
> # snp effect of solutions
> snp_eff = sol[2:4]
> snp_eff
[1] 0.6068652 -4.0802744 1.9341548
>
> # polygenic effect of solutions
> polge_eff = sol[5:30]
> polge_eff
[1] -6.311179e-15 1.220146e-01 -6.787405e-15 1.945577e-01 -1.042586e-01 9.507380e-02 -1.720163e-15 -1.656048e-15 -2.644430e-01 1.180100e-16
[11] -2.138676e-15 -6.615097e-16 -2.988163e-01 2.558718e-01 1.424285e-01 2.542361e-01 -8.517367e-02 2.705466e-01 -9.223368e-02 -1.807806e-01
[21] -7.584268e-15 1.279359e-01 1.279359e-01 1.279359e-01 1.279359e-01 1.279359e-01
>
> # 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,] 2.83345926
[2,] 0.29243931
[3,] 0.08047441
[4,] 3.55388367
[5,] -2.46054553
[6,] -3.78783512
[7,] 1.61972891
[8,] -2.46054553
[9,] 0.89930449
[10,] -0.31442586
[11,] -2.46054553
[12,] 0.29243931
[13,] -0.31442586
[14,] 2.22659408
>
> # GEBV : animal 13 ~ 26
> gebv = dgv + polge_eff[13:26]
> gebv
[,1]
[1,] 2.5346429
[2,] 0.5483112
[3,] 0.2229029
[4,] 3.8081198
[5,] -2.5457192
[6,] -3.5172885
[7,] 1.5274952
[8,] -2.6413261
[9,] 0.8993045
[10,] -0.1864899
[11,] -2.3326096
[12,] 0.4203752
[13,] -0.1864899
[14,] 2.3545300
>
SNP effect를 고정효과로 하여 추정한 후 개체의 SNP effect의 design matrix를 이용하여 DGV를 추정하였다. DGV와 polygenic effect를 합하여 GEBV를 계산하였다.
'Animal Breeding > R for Genetic Evaluation' 카테고리의 다른 글
| Mixed Linear Model(SNP-BLUP Model) for Computing SNP Effects, unweighted analysis (0) | 2020.12.19 |
|---|---|
| Fixed Effect Model for SNP Effects, weighted analysis (0) | 2020.12.16 |
| 임의 회귀 모형의 해를 이용한 개체 육종가 그래프 및 305일령 육종가 계산 (0) | 2020.11.26 |
| 다형질 모형의 차원 줄이기(Reduce the Dimension of Multivariate Models) (0) | 2020.09.29 |
| R을 이용한 일반적인 모형의 육종가 구하기 프로그래밍 (0) | 2020.08.25 |