Let's Study Animal Breeding

전체 글

R을 이용한 사회적 관계 모형(Social Interaction Model)의 육종가 구하기 2021.04.15
H 행렬, H 행렬의 역행렬 유도 2021.01.19
blupf90으로 single-step GBLUP Model 풀기 2021.01.05
Single-step GBLUP(ssGBLUP) 2021.01.05
GBLUP Model with residual polygenic effects 2020.12.29

R을 이용한 사회적 관계 모형(Social Interaction Model)의 육종가 구하기

투정이 2021. 4. 15. 17:20

2021. 4. 15. 17:20

# Social Interaction Model

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

# Raphael Mrode

# Example 8.1 p125

하나의 펜(pen) 또는 케이지에 여러 마리를 검정할 경우 개체들 사이의 경쟁과 협력에 따라서 검정 결과가 달라진다.

기본적으로 주어지는 자료는 다음과 같다.

7 1 4 1 1 5.5
8 1 4 1 2 9.8
9 2 5 1 2 4.9
10 1 4 2 1 8.23
11 2 5 2 2 7.5
12 3 6 2 2 10
13 2 5 3 1 4.5
14 3 6 3 2 8.4
15 3 6 3 1 6.4

1열 : 개체

2열 : 아비

3열 : 어미

4열 : 펜

5열 : 성별

6열 : 성장율

그런데 같은 펜에 있는 동기축에 대한 고려가 있어야 하므로 같은 펜 내의 개체를 제외한 다른 개체를 열로 추가하여야 한다. 동기축을 열에 추가한 결과는 다음과 같다.

7 1 4 1 1 5.5 8 9
8 1 4 1 2 9.8 7 9
9 2 5 1 2 4.9 7 8
10 1 4 2 1 8.23 11 12
11 2 5 2 2 7.5 10 12
12 3 6 2 2 10 10 11
13 2 5 3 1 4.5 14 15
14 3 6 3 2 8.4 13 15
15 3 6 3 1 6.4 13 14

7열 : 첫째 동기축(경쟁자)

8열 : 둘째 동기축(경쟁자)

모형

1. 성별효과 : 고정효과

2. 직접 육종가(direct EBV)

3. 조합 육종가(associative EBV) - 표현형을 동기축에 연결

4. 그룹(pen) 효과 : 임의 효과

5. 공통 환경 효과 - full-sib 효과 : 임의 효과

등을 포함한다.

모수

직접 육종가와 조합육종가의 분산 : 25.7, 3.6

직접 육종가와 조합육종가 사이의 공분산 : 2.25

공통 환경 분산 : 12.5

직접 효과의 잔차 분산 : 40.6

조합 효과의 잔차 분산 : 10.0

같은 펜에 있는 개체 사이의 상관 : 0.2

혈통

R 프로그램

# Social interaction effect model - Date : 2021-4-14

# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition
# Raphael Mrode
# Example 8.1 p125

# 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()

# pedigree and data

# read pedigree : animal sire dam
pedi = read.table("social_pedi.txt", header = FALSE, sep = "", col.names = c("animal", "sire", "dam"),stringsAsFactors = FALSE)
pedi = pedi[order(pedi$animal),]

# print
pedi

# read data : animal sire dam pen sex growth c1 c2
data = read.table("social_data.txt", header = FALSE, sep = "", col.names = c("animal", "sire", "dam", "pen", "sex", "growth", "c1", "c2"), stringsAsFactors = FALSE)

data

# growth, testing result
y = data[, 'growth']
y

# number of animals in the pen
nap = 3
nap

# genetic variances for direct and associative, covariance between them
G = matrix(c(25.7, 2.25, 2.25, 3.6), 2, 2)
G

# common environmental variance
sigma_c = 12.5
sigma_c

# residual variance variance for direct and associative
sigma_ed = 40.6
sigma_es = 10
sigma_ed
sigma_es

# correlation among pigs in the same pen(rho)
rho = 0.2
rho

# original residual variance
sigma_e = sigma_ed + (nap - 1) * sigma_es
sigma_e

# group variance
sigma_g = rho * sigma_e
sigma_g

# final residual variance
sigma_ef = sigma_e - sigma_g
sigma_ef

# alpha : variance component ratio of G
alpha = ginv(G) * sigma_ef
alpha

alpha1 = alpha[1,1]
alpha1

alpha2 = alpha[1,2]
alpha2

alpha3 = alpha[2,2]
alpha3

alpha4 = sigma_ef / sigma_g
alpha4

alpha5 = sigma_ef / sigma_c
alpha5

# design matrix : fixed effect
X = desgn(data[,'sex'])
X

# design matrix : direct effect
Zd = desgn(data[, 'animal'])
Zd

# design matrix : associative effect
others = data[,c("c1","c2")]
others

Zs = matrix(0, 9, 15)

for (i in 1:9) {
    Zs[i, others[i, 1]] = 1
    Zs[i, others[i, 2]] = 1
}

Zs

# design matrix : group(pen) effect
V = desgn(data[, 'pen'])
V

# number of group
no_g = ncol(V)
no_g

# design matrix : full-sib common environmental effect
ce = paste0(data[, 'sire'], '_', data[, 'dam'])
ce

W = desgn(ce)
W

# number of common environmental effects
no_ce = ncol(W)
no_ce
  

# Inverse matrix of Numerator Relationship Matrix
ai = ainv(pedi)
ai

# Left hand side
XPX  = t(X) %*% X
XPZd = t(X) %*% Zd
XPZs = t(X) %*% Zs
XPV  = t(X) %*% V
XPW  = t(X) %*% W

ZdPZd = t(Zd) %*% Zd
ZdPZs = t(Zd) %*% Zs
ZdPV  = t(Zd) %*% V
ZdPW  = t(Zd) %*% W

ZsPZs = t(Zs) %*% Zs
ZsPV  = t(Zs) %*% V
ZsPW  = t(Zs) %*% W

VPV  = t(V) %*% V
VPW  = t(V) %*% W

WPW  = t(W) %*% W

XPX
XPZd
XPZs
XPV
XPW

ZdPZd
ZdPZs
ZdPV
ZdPW

ZsPZs
ZsPV
ZsPW

VPV
VPW

WPW

LHS = rbind(
        cbind(XPX,     XPZd,                   XPZs,                  XPV,                       XPW),
        cbind(t(XPZd), ZdPZd    + ai * alpha1, ZdPZs + ai * alpha2,   ZdPV,                      ZdPW),
        cbind(t(XPZs), t(ZdPZs) + ai * alpha2, ZsPZs + ai * alpha3,   ZsPV,                      ZsPW),
        cbind(t(XPV),  t(ZdPV),                t(ZsPV),               VPV + diag(no_g) * alpha4, VPW),
        cbind(t(XPW),  t(ZdPW),                t(ZsPW),               t(VPW),                    WPW + diag(no_ce) * alpha5)
)

round(LHS,2)

# Right hand side
XPy  = t(X)  %*% y
ZdPy = t(Zd) %*% y
ZsPy = t(Zs) %*% y
VPy  = t(V)  %*% y
WPy  = t(W)  %*% y

XPy
ZdPy
ZsPy
VPy
WPy

RHS = rbind(XPy,
            ZdPy,
            ZsPy,
            VPy,
            WPy
            )

RHS

# generalized inverse of LHS
gi_LHS = ginv(LHS)
round(gi_LHS, 3)

diag(gi_LHS)

# solution
sol = gi_LHS %*% RHS

# fixed effect : sex effect
sol_f1 = sol[1:2]
sol_f1

# direct EBV
sol_dbv = sol[3:17]

# associative EBV
sol_sbv = sol [18:32]

sol_bv = cbind(sol_dbv, sol_sbv, sol_dbv + sol_sbv)
sol_bv

# random effect : group effect
sol_r1 = sol[33:35]
sol_r1

# random effect : common environmental effect
sol_r2 = sol[36:38]
sol_r2

RStudio 실행 화면

실행 결과

> # Social interaction effect model - Date : 2021-4-14
> 
> # Linear Models for the Prediction of Animal Breeding Values, 3rd Edition
> # Raphael Mrode
> # Example 8.1 p125
> 
> # 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/2021/job/공부하기/01_Linear Models for the Prediction of Animal Breeding Values_3rd_Mrode/28_social interaction model_R"
> 
> # pedigree and data
> 
> # read pedigree : animal sire dam
> pedi = read.table("social_pedi.txt", header = FALSE, sep = "", col.names = c("animal", "sire", "dam"),stringsAsFactors = FALSE)
> pedi = pedi[order(pedi$animal),]
> 
> # print
> pedi
   animal sire dam
1       1    0   0
2       2    0   0
3       3    0   0
4       4    0   0
5       5    0   0
6       6    0   0
7       7    1   4
8       8    1   4
9       9    2   5
10     10    1   4
11     11    2   5
12     12    3   6
13     13    2   5
14     14    3   6
15     15    3   6
> 
> # read data : animal sire dam pen sex growth c1 c2
> data = read.table("social_data.txt", header = FALSE, sep = "", col.names = c("animal", "sire", "dam", "pen", "sex", "growth", "c1", "c2"), stringsAsFactors = FALSE)
> 
> data
  animal sire dam pen sex growth c1 c2
1      7    1   4   1   1   5.50  8  9
2      8    1   4   1   2   9.80  7  9
3      9    2   5   1   2   4.90  7  8
4     10    1   4   2   1   8.23 11 12
5     11    2   5   2   2   7.50 10 12
6     12    3   6   2   2  10.00 10 11
7     13    2   5   3   1   4.50 14 15
8     14    3   6   3   2   8.40 13 15
9     15    3   6   3   1   6.40 13 14
> 
> # growth, testing result
> y = data[, 'growth']
> y
[1]  5.50  9.80  4.90  8.23  7.50 10.00  4.50  8.40  6.40
> 
> # number of animals in the pen
> nap = 3
> nap
[1] 3
> 
> # genetic variances for direct and associative, covariance between them
> G = matrix(c(25.7, 2.25, 2.25, 3.6), 2, 2)
> G
      [,1] [,2]
[1,] 25.70 2.25
[2,]  2.25 3.60
> 
> # common environmental variance
> sigma_c = 12.5
> sigma_c
[1] 12.5
> 
> # residual variance variance for direct and associative
> sigma_ed = 40.6
> sigma_es = 10
> sigma_ed
[1] 40.6
> sigma_es
[1] 10
> 
> # correlation among pigs in the same pen(rho)
> rho = 0.2
> rho
[1] 0.2
> 
> # original residual variance
> sigma_e = sigma_ed + (nap - 1) * sigma_es
> sigma_e
[1] 60.6
> 
> # group variance
> sigma_g = rho * sigma_e
> sigma_g
[1] 12.12
> 
> # final residual variance
> sigma_ef = sigma_e - sigma_g
> sigma_ef
[1] 48.48
> 
> # alpha : variance component ratio of G
> alpha = ginv(G) * sigma_ef
> alpha
          [,1]      [,2]
[1,]  1.995575 -1.247234
[2,] -1.247234 14.246188
> 
> alpha1 = alpha[1,1]
> alpha1
[1] 1.995575
> 
> alpha2 = alpha[1,2]
> alpha2
[1] -1.247234
> 
> alpha3 = alpha[2,2]
> alpha3
[1] 14.24619
> 
> alpha4 = sigma_ef / sigma_g
> alpha4
[1] 4
> 
> alpha5 = sigma_ef / sigma_c
> alpha5
[1] 3.8784
> 
> # design matrix : fixed effect
> X = desgn(data[,'sex'])
> X
      [,1] [,2]
 [1,]    1    0
 [2,]    0    1
 [3,]    0    1
 [4,]    1    0
 [5,]    0    1
 [6,]    0    1
 [7,]    1    0
 [8,]    0    1
 [9,]    1    0
> 
> # design matrix : direct effect
> Zd = desgn(data[, 'animal'])
> Zd
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
 [1,]    0    0    0    0    0    0    1    0    0     0     0     0     0     0     0
 [2,]    0    0    0    0    0    0    0    1    0     0     0     0     0     0     0
 [3,]    0    0    0    0    0    0    0    0    1     0     0     0     0     0     0
 [4,]    0    0    0    0    0    0    0    0    0     1     0     0     0     0     0
 [5,]    0    0    0    0    0    0    0    0    0     0     1     0     0     0     0
 [6,]    0    0    0    0    0    0    0    0    0     0     0     1     0     0     0
 [7,]    0    0    0    0    0    0    0    0    0     0     0     0     1     0     0
 [8,]    0    0    0    0    0    0    0    0    0     0     0     0     0     1     0
 [9,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     1
> 
> # design matrix : associative effect
> others = data[,c("c1","c2")]
> others
  c1 c2
1  8  9
2  7  9
3  7  8
4 11 12
5 10 12
6 10 11
7 14 15
8 13 15
9 13 14
> 
> Zs = matrix(0, 9, 15)
> 
> for (i in 1:9) {
+     Zs[i, others[i, 1]] = 1
+     Zs[i, others[i, 2]] = 1
+ }
> 
> Zs
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
 [1,]    0    0    0    0    0    0    0    1    1     0     0     0     0     0     0
 [2,]    0    0    0    0    0    0    1    0    1     0     0     0     0     0     0
 [3,]    0    0    0    0    0    0    1    1    0     0     0     0     0     0     0
 [4,]    0    0    0    0    0    0    0    0    0     0     1     1     0     0     0
 [5,]    0    0    0    0    0    0    0    0    0     1     0     1     0     0     0
 [6,]    0    0    0    0    0    0    0    0    0     1     1     0     0     0     0
 [7,]    0    0    0    0    0    0    0    0    0     0     0     0     0     1     1
 [8,]    0    0    0    0    0    0    0    0    0     0     0     0     1     0     1
 [9,]    0    0    0    0    0    0    0    0    0     0     0     0     1     1     0
> 
> # design matrix : group(pen) effect
> V = desgn(data[, 'pen'])
> V
      [,1] [,2] [,3]
 [1,]    1    0    0
 [2,]    1    0    0
 [3,]    1    0    0
 [4,]    0    1    0
 [5,]    0    1    0
 [6,]    0    1    0
 [7,]    0    0    1
 [8,]    0    0    1
 [9,]    0    0    1
> 
> # number of group
> no_g = ncol(V)
> no_g
[1] 3
> 
> # design matrix : full-sib common environmental effect
> ce = paste0(data[, 'sire'], '_', data[, 'dam'])
> ce
[1] "1_4" "1_4" "2_5" "1_4" "2_5" "3_6" "2_5" "3_6" "3_6"
> 
> W = desgn(ce)
> W
      [,1] [,2] [,3]
 [1,]    1    0    0
 [2,]    1    0    0
 [3,]    0    1    0
 [4,]    1    0    0
 [5,]    0    1    0
 [6,]    0    0    1
 [7,]    0    1    0
 [8,]    0    0    1
 [9,]    0    0    1
> 
> # number of common environmental effects
> no_ce = ncol(W)
> no_ce
[1] 3
>   
> 
> # Inverse matrix of Numerator Relationship Matrix
> ai = ainv(pedi)
> ai
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
 [1,]  2.5  0.0  0.0  1.5  0.0  0.0   -1   -1    0    -1     0     0     0     0     0
 [2,]  0.0  2.5  0.0  0.0  1.5  0.0    0    0   -1     0    -1     0    -1     0     0
 [3,]  0.0  0.0  2.5  0.0  0.0  1.5    0    0    0     0     0    -1     0    -1    -1
 [4,]  1.5  0.0  0.0  2.5  0.0  0.0   -1   -1    0    -1     0     0     0     0     0
 [5,]  0.0  1.5  0.0  0.0  2.5  0.0    0    0   -1     0    -1     0    -1     0     0
 [6,]  0.0  0.0  1.5  0.0  0.0  2.5    0    0    0     0     0    -1     0    -1    -1
 [7,] -1.0  0.0  0.0 -1.0  0.0  0.0    2    0    0     0     0     0     0     0     0
 [8,] -1.0  0.0  0.0 -1.0  0.0  0.0    0    2    0     0     0     0     0     0     0
 [9,]  0.0 -1.0  0.0  0.0 -1.0  0.0    0    0    2     0     0     0     0     0     0
[10,] -1.0  0.0  0.0 -1.0  0.0  0.0    0    0    0     2     0     0     0     0     0
[11,]  0.0 -1.0  0.0  0.0 -1.0  0.0    0    0    0     0     2     0     0     0     0
[12,]  0.0  0.0 -1.0  0.0  0.0 -1.0    0    0    0     0     0     2     0     0     0
[13,]  0.0 -1.0  0.0  0.0 -1.0  0.0    0    0    0     0     0     0     2     0     0
[14,]  0.0  0.0 -1.0  0.0  0.0 -1.0    0    0    0     0     0     0     0     2     0
[15,]  0.0  0.0 -1.0  0.0  0.0 -1.0    0    0    0     0     0     0     0     0     2
> 
> # Left hand side
> XPX  = t(X) %*% X
> XPZd = t(X) %*% Zd
> XPZs = t(X) %*% Zs
> XPV  = t(X) %*% V
> XPW  = t(X) %*% W
> 
> ZdPZd = t(Zd) %*% Zd
> ZdPZs = t(Zd) %*% Zs
> ZdPV  = t(Zd) %*% V
> ZdPW  = t(Zd) %*% W
> 
> ZsPZs = t(Zs) %*% Zs
> ZsPV  = t(Zs) %*% V
> ZsPW  = t(Zs) %*% W
> 
> VPV  = t(V) %*% V
> VPW  = t(V) %*% W
> 
> WPW  = t(W) %*% W
> 
> XPX
     [,1] [,2]
[1,]    4    0
[2,]    0    5
> XPZd
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
[1,]    0    0    0    0    0    0    1    0    0     1     0     0     1     0     1
[2,]    0    0    0    0    0    0    0    1    1     0     1     1     0     1     0
> XPZs
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
[1,]    0    0    0    0    0    0    0    1    1     0     1     1     1     2     1
[2,]    0    0    0    0    0    0    2    1    1     2     1     1     1     0     1
> XPV
     [,1] [,2] [,3]
[1,]    1    1    2
[2,]    2    2    1
> XPW
     [,1] [,2] [,3]
[1,]    2    1    1
[2,]    1    2    2
> 
> ZdPZd
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
 [1,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0
 [2,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0
 [3,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0
 [4,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0
 [5,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0
 [6,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0
 [7,]    0    0    0    0    0    0    1    0    0     0     0     0     0     0     0
 [8,]    0    0    0    0    0    0    0    1    0     0     0     0     0     0     0
 [9,]    0    0    0    0    0    0    0    0    1     0     0     0     0     0     0
[10,]    0    0    0    0    0    0    0    0    0     1     0     0     0     0     0
[11,]    0    0    0    0    0    0    0    0    0     0     1     0     0     0     0
[12,]    0    0    0    0    0    0    0    0    0     0     0     1     0     0     0
[13,]    0    0    0    0    0    0    0    0    0     0     0     0     1     0     0
[14,]    0    0    0    0    0    0    0    0    0     0     0     0     0     1     0
[15,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     1
> ZdPZs
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
 [1,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0
 [2,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0
 [3,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0
 [4,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0
 [5,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0
 [6,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0
 [7,]    0    0    0    0    0    0    0    1    1     0     0     0     0     0     0
 [8,]    0    0    0    0    0    0    1    0    1     0     0     0     0     0     0
 [9,]    0    0    0    0    0    0    1    1    0     0     0     0     0     0     0
[10,]    0    0    0    0    0    0    0    0    0     0     1     1     0     0     0
[11,]    0    0    0    0    0    0    0    0    0     1     0     1     0     0     0
[12,]    0    0    0    0    0    0    0    0    0     1     1     0     0     0     0
[13,]    0    0    0    0    0    0    0    0    0     0     0     0     0     1     1
[14,]    0    0    0    0    0    0    0    0    0     0     0     0     1     0     1
[15,]    0    0    0    0    0    0    0    0    0     0     0     0     1     1     0
> ZdPV
      [,1] [,2] [,3]
 [1,]    0    0    0
 [2,]    0    0    0
 [3,]    0    0    0
 [4,]    0    0    0
 [5,]    0    0    0
 [6,]    0    0    0
 [7,]    1    0    0
 [8,]    1    0    0
 [9,]    1    0    0
[10,]    0    1    0
[11,]    0    1    0
[12,]    0    1    0
[13,]    0    0    1
[14,]    0    0    1
[15,]    0    0    1
> ZdPW
      [,1] [,2] [,3]
 [1,]    0    0    0
 [2,]    0    0    0
 [3,]    0    0    0
 [4,]    0    0    0
 [5,]    0    0    0
 [6,]    0    0    0
 [7,]    1    0    0
 [8,]    1    0    0
 [9,]    0    1    0
[10,]    1    0    0
[11,]    0    1    0
[12,]    0    0    1
[13,]    0    1    0
[14,]    0    0    1
[15,]    0    0    1
> 
> ZsPZs
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
 [1,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0
 [2,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0
 [3,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0
 [4,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0
 [5,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0
 [6,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0
 [7,]    0    0    0    0    0    0    2    1    1     0     0     0     0     0     0
 [8,]    0    0    0    0    0    0    1    2    1     0     0     0     0     0     0
 [9,]    0    0    0    0    0    0    1    1    2     0     0     0     0     0     0
[10,]    0    0    0    0    0    0    0    0    0     2     1     1     0     0     0
[11,]    0    0    0    0    0    0    0    0    0     1     2     1     0     0     0
[12,]    0    0    0    0    0    0    0    0    0     1     1     2     0     0     0
[13,]    0    0    0    0    0    0    0    0    0     0     0     0     2     1     1
[14,]    0    0    0    0    0    0    0    0    0     0     0     0     1     2     1
[15,]    0    0    0    0    0    0    0    0    0     0     0     0     1     1     2
> ZsPV
      [,1] [,2] [,3]
 [1,]    0    0    0
 [2,]    0    0    0
 [3,]    0    0    0
 [4,]    0    0    0
 [5,]    0    0    0
 [6,]    0    0    0
 [7,]    2    0    0
 [8,]    2    0    0
 [9,]    2    0    0
[10,]    0    2    0
[11,]    0    2    0
[12,]    0    2    0
[13,]    0    0    2
[14,]    0    0    2
[15,]    0    0    2
> ZsPW
      [,1] [,2] [,3]
 [1,]    0    0    0
 [2,]    0    0    0
 [3,]    0    0    0
 [4,]    0    0    0
 [5,]    0    0    0
 [6,]    0    0    0
 [7,]    1    1    0
 [8,]    1    1    0
 [9,]    2    0    0
[10,]    0    1    1
[11,]    1    0    1
[12,]    1    1    0
[13,]    0    0    2
[14,]    0    1    1
[15,]    0    1    1
> 
> VPV
     [,1] [,2] [,3]
[1,]    3    0    0
[2,]    0    3    0
[3,]    0    0    3
> VPW
     [,1] [,2] [,3]
[1,]    2    1    0
[2,]    1    1    1
[3,]    0    1    2
> 
> WPW
     [,1] [,2] [,3]
[1,]    3    0    0
[2,]    0    3    0
[3,]    0    0    3
> 
> LHS = rbind(
+         cbind(XPX,     XPZd,                   XPZs,                  XPV,                       XPW),
+         cbind(t(XPZd), ZdPZd    + ai * alpha1, ZdPZs + ai * alpha2,   ZdPV,                      ZdPW),
+         cbind(t(XPZs), t(ZdPZs) + ai * alpha2, ZsPZs + ai * alpha3,   ZsPV,                      ZsPW),
+         cbind(t(XPV),  t(ZdPV),                t(ZsPV),               VPV + diag(no_g) * alpha4, VPW),
+         cbind(t(XPW),  t(ZdPW),                t(ZsPW),               t(VPW),                    WPW + diag(no_ce) * alpha5)
+ )
> 
> round(LHS,2)
      [,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]  [,28]  [,29]
 [1,]    4    0  0.00  0.00  0.00  0.00  0.00  0.00  1.00  0.00  0.00  1.00  0.00  0.00  1.00  0.00  1.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   1.00   1.00   0.00   1.00   1.00
 [2,]    0    5  0.00  0.00  0.00  0.00  0.00  0.00  0.00  1.00  1.00  0.00  1.00  1.00  0.00  1.00  0.00   0.00   0.00   0.00   0.00   0.00   0.00   2.00   1.00   1.00   2.00   1.00   1.00
 [3,]    0    0  4.99  0.00  0.00  2.99  0.00  0.00 -2.00 -2.00  0.00 -2.00  0.00  0.00  0.00  0.00  0.00  -3.12   0.00   0.00  -1.87   0.00   0.00   1.25   1.25   0.00   1.25   0.00   0.00
 [4,]    0    0  0.00  4.99  0.00  0.00  2.99  0.00  0.00  0.00 -2.00  0.00 -2.00  0.00 -2.00  0.00  0.00   0.00  -3.12   0.00   0.00  -1.87   0.00   0.00   0.00   1.25   0.00   1.25   0.00
 [5,]    0    0  0.00  0.00  4.99  0.00  0.00  2.99  0.00  0.00  0.00  0.00  0.00 -2.00  0.00 -2.00 -2.00   0.00   0.00  -3.12   0.00   0.00  -1.87   0.00   0.00   0.00   0.00   0.00   1.25
 [6,]    0    0  2.99  0.00  0.00  4.99  0.00  0.00 -2.00 -2.00  0.00 -2.00  0.00  0.00  0.00  0.00  0.00  -1.87   0.00   0.00  -3.12   0.00   0.00   1.25   1.25   0.00   1.25   0.00   0.00
 [7,]    0    0  0.00  2.99  0.00  0.00  4.99  0.00  0.00  0.00 -2.00  0.00 -2.00  0.00 -2.00  0.00  0.00   0.00  -1.87   0.00   0.00  -3.12   0.00   0.00   0.00   1.25   0.00   1.25   0.00
 [8,]    0    0  0.00  0.00  2.99  0.00  0.00  4.99  0.00  0.00  0.00  0.00  0.00 -2.00  0.00 -2.00 -2.00   0.00   0.00  -1.87   0.00   0.00  -3.12   0.00   0.00   0.00   0.00   0.00   1.25
 [9,]    1    0 -2.00  0.00  0.00 -2.00  0.00  0.00  4.99  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00   1.25   0.00   0.00   1.25   0.00   0.00  -2.49   1.00   1.00   0.00   0.00   0.00
[10,]    0    1 -2.00  0.00  0.00 -2.00  0.00  0.00  0.00  4.99  0.00  0.00  0.00  0.00  0.00  0.00  0.00   1.25   0.00   0.00   1.25   0.00   0.00   1.00  -2.49   1.00   0.00   0.00   0.00
[11,]    0    1  0.00 -2.00  0.00  0.00 -2.00  0.00  0.00  0.00  4.99  0.00  0.00  0.00  0.00  0.00  0.00   0.00   1.25   0.00   0.00   1.25   0.00   1.00   1.00  -2.49   0.00   0.00   0.00
[12,]    1    0 -2.00  0.00  0.00 -2.00  0.00  0.00  0.00  0.00  0.00  4.99  0.00  0.00  0.00  0.00  0.00   1.25   0.00   0.00   1.25   0.00   0.00   0.00   0.00   0.00  -2.49   1.00   1.00
[13,]    0    1  0.00 -2.00  0.00  0.00 -2.00  0.00  0.00  0.00  0.00  0.00  4.99  0.00  0.00  0.00  0.00   0.00   1.25   0.00   0.00   1.25   0.00   0.00   0.00   0.00   1.00  -2.49   1.00
[14,]    0    1  0.00  0.00 -2.00  0.00  0.00 -2.00  0.00  0.00  0.00  0.00  0.00  4.99  0.00  0.00  0.00   0.00   0.00   1.25   0.00   0.00   1.25   0.00   0.00   0.00   1.00   1.00  -2.49
[15,]    1    0  0.00 -2.00  0.00  0.00 -2.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  4.99  0.00  0.00   0.00   1.25   0.00   0.00   1.25   0.00   0.00   0.00   0.00   0.00   0.00   0.00
[16,]    0    1  0.00  0.00 -2.00  0.00  0.00 -2.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  4.99  0.00   0.00   0.00   1.25   0.00   0.00   1.25   0.00   0.00   0.00   0.00   0.00   0.00
[17,]    1    0  0.00  0.00 -2.00  0.00  0.00 -2.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  0.00  4.99   0.00   0.00   1.25   0.00   0.00   1.25   0.00   0.00   0.00   0.00   0.00   0.00
[18,]    0    0 -3.12  0.00  0.00 -1.87  0.00  0.00  1.25  1.25  0.00  1.25  0.00  0.00  0.00  0.00  0.00  35.62   0.00   0.00  21.37   0.00   0.00 -14.25 -14.25   0.00 -14.25   0.00   0.00
[19,]    0    0  0.00 -3.12  0.00  0.00 -1.87  0.00  0.00  0.00  1.25  0.00  1.25  0.00  1.25  0.00  0.00   0.00  35.62   0.00   0.00  21.37   0.00   0.00   0.00 -14.25   0.00 -14.25   0.00
[20,]    0    0  0.00  0.00 -3.12  0.00  0.00 -1.87  0.00  0.00  0.00  0.00  0.00  1.25  0.00  1.25  1.25   0.00   0.00  35.62   0.00   0.00  21.37   0.00   0.00   0.00   0.00   0.00 -14.25
[21,]    0    0 -1.87  0.00  0.00 -3.12  0.00  0.00  1.25  1.25  0.00  1.25  0.00  0.00  0.00  0.00  0.00  21.37   0.00   0.00  35.62   0.00   0.00 -14.25 -14.25   0.00 -14.25   0.00   0.00
[22,]    0    0  0.00 -1.87  0.00  0.00 -3.12  0.00  0.00  0.00  1.25  0.00  1.25  0.00  1.25  0.00  0.00   0.00  21.37   0.00   0.00  35.62   0.00   0.00   0.00 -14.25   0.00 -14.25   0.00
[23,]    0    0  0.00  0.00 -1.87  0.00  0.00 -3.12  0.00  0.00  0.00  0.00  0.00  1.25  0.00  1.25  1.25   0.00   0.00  21.37   0.00   0.00  35.62   0.00   0.00   0.00   0.00   0.00 -14.25
[24,]    0    2  1.25  0.00  0.00  1.25  0.00  0.00 -2.49  1.00  1.00  0.00  0.00  0.00  0.00  0.00  0.00 -14.25   0.00   0.00 -14.25   0.00   0.00  30.49   1.00   1.00   0.00   0.00   0.00
[25,]    1    1  1.25  0.00  0.00  1.25  0.00  0.00  1.00 -2.49  1.00  0.00  0.00  0.00  0.00  0.00  0.00 -14.25   0.00   0.00 -14.25   0.00   0.00   1.00  30.49   1.00   0.00   0.00   0.00
[26,]    1    1  0.00  1.25  0.00  0.00  1.25  0.00  1.00  1.00 -2.49  0.00  0.00  0.00  0.00  0.00  0.00   0.00 -14.25   0.00   0.00 -14.25   0.00   1.00   1.00  30.49   0.00   0.00   0.00
       [,30]  [,31]  [,32] [,33] [,34] [,35] [,36] [,37] [,38]
 [1,]   1.00   2.00   1.00     1     1     2  2.00  1.00  1.00
 [2,]   1.00   0.00   1.00     2     2     1  1.00  2.00  2.00
 [3,]   0.00   0.00   0.00     0     0     0  0.00  0.00  0.00
 [4,]   1.25   0.00   0.00     0     0     0  0.00  0.00  0.00
 [5,]   0.00   1.25   1.25     0     0     0  0.00  0.00  0.00
 [6,]   0.00   0.00   0.00     0     0     0  0.00  0.00  0.00
 [7,]   1.25   0.00   0.00     0     0     0  0.00  0.00  0.00
 [8,]   0.00   1.25   1.25     0     0     0  0.00  0.00  0.00
 [9,]   0.00   0.00   0.00     1     0     0  1.00  0.00  0.00
[10,]   0.00   0.00   0.00     1     0     0  1.00  0.00  0.00
[11,]   0.00   0.00   0.00     1     0     0  0.00  1.00  0.00
[12,]   0.00   0.00   0.00     0     1     0  1.00  0.00  0.00
[13,]   0.00   0.00   0.00     0     1     0  0.00  1.00  0.00
[14,]   0.00   0.00   0.00     0     1     0  0.00  0.00  1.00
[15,]  -2.49   1.00   1.00     0     0     1  0.00  1.00  0.00
[16,]   1.00  -2.49   1.00     0     0     1  0.00  0.00  1.00
[17,]   1.00   1.00  -2.49     0     0     1  0.00  0.00  1.00
[18,]   0.00   0.00   0.00     0     0     0  0.00  0.00  0.00
[19,] -14.25   0.00   0.00     0     0     0  0.00  0.00  0.00
[20,]   0.00 -14.25 -14.25     0     0     0  0.00  0.00  0.00
[21,]   0.00   0.00   0.00     0     0     0  0.00  0.00  0.00
[22,] -14.25   0.00   0.00     0     0     0  0.00  0.00  0.00
[23,]   0.00 -14.25 -14.25     0     0     0  0.00  0.00  0.00
[24,]   0.00   0.00   0.00     2     0     0  1.00  1.00  0.00
[25,]   0.00   0.00   0.00     2     0     0  1.00  1.00  0.00
[26,]   0.00   0.00   0.00     2     0     0  2.00  0.00  0.00
 [ getOption("max.print") 에 도달했습니다 -- 12 행들을 생략합니다 ]
> 
> # Right hand side
> XPy  = t(X)  %*% y
> ZdPy = t(Zd) %*% y
> ZsPy = t(Zs) %*% y
> VPy  = t(V)  %*% y
> WPy  = t(W)  %*% y
> 
> XPy
      [,1]
[1,] 24.63
[2,] 40.60
> ZdPy
       [,1]
 [1,]  0.00
 [2,]  0.00
 [3,]  0.00
 [4,]  0.00
 [5,]  0.00
 [6,]  0.00
 [7,]  5.50
 [8,]  9.80
 [9,]  4.90
[10,]  8.23
[11,]  7.50
[12,] 10.00
[13,]  4.50
[14,]  8.40
[15,]  6.40
> ZsPy
       [,1]
 [1,]  0.00
 [2,]  0.00
 [3,]  0.00
 [4,]  0.00
 [5,]  0.00
 [6,]  0.00
 [7,] 14.70
 [8,] 10.40
 [9,] 15.30
[10,] 17.50
[11,] 18.23
[12,] 15.73
[13,] 14.80
[14,] 10.90
[15,] 12.90
> VPy
      [,1]
[1,] 20.20
[2,] 25.73
[3,] 19.30
> WPy
      [,1]
[1,] 23.53
[2,] 16.90
[3,] 24.80
> 
> RHS = rbind(XPy,
+             ZdPy,
+             ZsPy,
+             VPy,
+             WPy
+             )
> 
> RHS
       [,1]
 [1,] 24.63
 [2,] 40.60
 [3,]  0.00
 [4,]  0.00
 [5,]  0.00
 [6,]  0.00
 [7,]  0.00
 [8,]  0.00
 [9,]  5.50
[10,]  9.80
[11,]  4.90
[12,]  8.23
[13,]  7.50
[14,] 10.00
[15,]  4.50
[16,]  8.40
[17,]  6.40
[18,]  0.00
[19,]  0.00
[20,]  0.00
[21,]  0.00
[22,]  0.00
[23,]  0.00
[24,] 14.70
[25,] 10.40
[26,] 15.30
[27,] 17.50
[28,] 18.23
[29,] 15.73
[30,] 14.80
[31,] 10.90
[32,] 12.90
[33,] 20.20
[34,] 25.73
[35,] 19.30
[36,] 23.53
[37,] 16.90
[38,] 24.80
> 
> # 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]
 [1,]  0.712  0.340 -0.129 -0.093 -0.089 -0.129 -0.093 -0.089 -0.197 -0.126 -0.107 -0.193 -0.102 -0.103 -0.163 -0.091 -0.163 -0.023 -0.033 -0.042 -0.023 -0.033 -0.042 -0.029 -0.033 -0.042
 [2,]  0.340  0.638 -0.084 -0.112 -0.115 -0.084 -0.112 -0.115 -0.091 -0.149 -0.164 -0.095 -0.167 -0.167 -0.119 -0.176 -0.119 -0.040 -0.032 -0.025 -0.040 -0.032 -0.025 -0.055 -0.052 -0.044
 [3,] -0.129 -0.084  0.489  0.020  0.022 -0.041  0.020  0.022  0.213  0.204  0.024  0.213  0.024  0.027  0.031  0.025  0.034  0.042 -0.001  0.005 -0.004 -0.001  0.005  0.016  0.017 -0.005
 [4,] -0.093 -0.112  0.020  0.489  0.022  0.020 -0.041  0.022  0.025  0.029  0.211  0.025  0.211  0.030  0.208  0.030  0.027 -0.001  0.049 -0.002 -0.001  0.002 -0.002 -0.001 -0.001  0.026
 [5,] -0.089 -0.115  0.022  0.022  0.487  0.022  0.022 -0.043  0.027  0.032  0.030  0.027  0.030  0.209  0.026  0.209  0.204  0.005 -0.001  0.043  0.005 -0.001 -0.004  0.008  0.007  0.001
 [6,] -0.129 -0.084 -0.041  0.020  0.022  0.489  0.020  0.022  0.213  0.204  0.024  0.213  0.024  0.027  0.031  0.025  0.034 -0.004 -0.001  0.005  0.042 -0.001  0.005  0.016  0.017 -0.005
 [7,] -0.093 -0.112  0.020 -0.041  0.022  0.020  0.489  0.022  0.025  0.029  0.211  0.025  0.211  0.030  0.208  0.030  0.027 -0.001  0.002 -0.002 -0.001  0.049 -0.002 -0.001 -0.001  0.026
 [8,] -0.089 -0.115  0.022  0.022 -0.043  0.022  0.022  0.487  0.027  0.032  0.030  0.027  0.030  0.209  0.026  0.209  0.204  0.005 -0.001 -0.004  0.005 -0.001  0.043  0.008  0.007  0.001
 [9,] -0.197 -0.091  0.213  0.025  0.027  0.213  0.025  0.027  0.438  0.199  0.029  0.216  0.027  0.031  0.044  0.028  0.049  0.015 -0.001  0.010  0.015 -0.001  0.010  0.034  0.009 -0.008
[10,] -0.126 -0.149  0.204  0.029  0.032  0.204  0.029  0.032  0.199  0.421  0.040  0.198  0.039  0.044  0.035  0.044  0.040  0.018 -0.002  0.007  0.018 -0.002  0.007  0.013  0.039 -0.008
[11,] -0.107 -0.164  0.024  0.211  0.030  0.024  0.211  0.030  0.029  0.040  0.432  0.027  0.211  0.044  0.200  0.044  0.033 -0.005  0.026  0.002 -0.005  0.026  0.002 -0.008 -0.009  0.049
[12,] -0.193 -0.095  0.213  0.025  0.027  0.213  0.025  0.027  0.216  0.198  0.027  0.437  0.029  0.032  0.045  0.029  0.048  0.020 -0.001  0.005  0.020 -0.001  0.005  0.018  0.020 -0.003
[13,] -0.102 -0.167  0.024  0.211  0.030  0.024  0.211  0.030  0.027  0.039  0.211  0.029  0.432  0.044  0.202  0.045  0.032  0.000  0.026 -0.003  0.000  0.026 -0.003  0.002  0.001  0.029
[14,] -0.103 -0.167  0.027  0.030  0.209  0.027  0.030  0.209  0.031  0.044  0.044  0.032  0.044  0.428  0.031  0.209  0.197  0.004 -0.002  0.021  0.004 -0.002  0.021  0.008  0.008  0.002
[15,] -0.163 -0.119  0.031  0.208  0.026  0.031  0.208  0.026  0.044  0.035  0.200  0.045  0.202  0.031  0.430  0.032  0.041  0.002  0.026 -0.005  0.002  0.026 -0.005  0.003  0.003  0.028
[16,] -0.091 -0.176  0.025  0.030  0.209  0.025  0.030  0.209  0.028  0.044  0.044  0.029  0.045  0.209  0.032  0.431  0.197  0.009 -0.002  0.015  0.009 -0.002  0.015  0.014  0.013  0.002
[17,] -0.163 -0.119  0.034  0.027  0.204  0.034  0.027  0.204  0.049  0.040  0.033  0.048  0.032  0.197  0.041  0.197  0.424  0.006 -0.002  0.019  0.006 -0.002  0.019  0.009  0.009  0.001
[18,] -0.023 -0.040  0.042 -0.001  0.005 -0.004 -0.001  0.005  0.015  0.018 -0.005  0.020  0.000  0.004  0.002  0.009  0.006  0.072  0.000  0.002 -0.002  0.000  0.002  0.034  0.034 -0.001
[19,] -0.033 -0.032 -0.001  0.049 -0.001 -0.001  0.002 -0.001 -0.001 -0.002  0.026 -0.001  0.026 -0.002  0.026 -0.002 -0.002  0.000  0.074  0.000  0.000  0.000  0.000  0.000  0.000  0.037
[20,] -0.042 -0.025  0.005 -0.002  0.043  0.005 -0.002 -0.004  0.010  0.007  0.002  0.005 -0.003  0.021 -0.005  0.015  0.019  0.002  0.000  0.072  0.002  0.000 -0.002  0.003  0.003  0.001
[21,] -0.023 -0.040 -0.004 -0.001  0.005  0.042 -0.001  0.005  0.015  0.018 -0.005  0.020  0.000  0.004  0.002  0.009  0.006 -0.002  0.000  0.002  0.072  0.000  0.002  0.034  0.034 -0.001
[22,] -0.033 -0.032 -0.001  0.002 -0.001 -0.001  0.049 -0.001 -0.001 -0.002  0.026 -0.001  0.026 -0.002  0.026 -0.002 -0.002  0.000  0.000  0.000  0.000  0.074  0.000  0.000  0.000  0.037
[23,] -0.042 -0.025  0.005 -0.002 -0.004  0.005 -0.002  0.043  0.010  0.007  0.002  0.005 -0.003  0.021 -0.005  0.015  0.019  0.002  0.000 -0.002  0.002  0.000  0.072  0.003  0.003  0.001
[24,] -0.029 -0.055  0.016 -0.001  0.008  0.016 -0.001  0.008  0.034  0.013 -0.008  0.018  0.002  0.008  0.003  0.014  0.009  0.034  0.000  0.003  0.034  0.000  0.003  0.070  0.032 -0.002
[25,] -0.033 -0.052  0.017 -0.001  0.007  0.017 -0.001  0.007  0.009  0.039 -0.009  0.020  0.001  0.008  0.003  0.013  0.009  0.034  0.000  0.003  0.034  0.000  0.003  0.032  0.069 -0.002
[26,] -0.042 -0.044 -0.005  0.026  0.001 -0.005  0.026  0.001 -0.008 -0.008  0.049 -0.003  0.029  0.002  0.028  0.002  0.001 -0.001  0.037  0.001 -0.001  0.037  0.001 -0.002 -0.002  0.073
       [,27]  [,28]  [,29]  [,30]  [,31]  [,32]  [,33]  [,34]  [,35]  [,36]  [,37]  [,38]
 [1,] -0.030 -0.043 -0.051 -0.047 -0.059 -0.055 -0.063 -0.070 -0.117 -0.118 -0.074 -0.065
 [2,] -0.054 -0.043 -0.037 -0.040 -0.030 -0.034 -0.099 -0.094 -0.056 -0.060 -0.095 -0.102
 [3,]  0.019 -0.002  0.005  0.002  0.009  0.008 -0.021 -0.002  0.022 -0.040  0.021  0.019
 [4,] -0.001  0.026 -0.002  0.026 -0.002 -0.002  0.001  0.001 -0.002  0.021 -0.044  0.023
 [5,]  0.005 -0.002  0.020 -0.005  0.017  0.017  0.020  0.001 -0.021  0.019  0.023 -0.042
 [6,]  0.019 -0.002  0.005  0.002  0.009  0.008 -0.021 -0.002  0.022 -0.040  0.021  0.019
 [7,] -0.001  0.026 -0.002  0.026 -0.002 -0.002  0.001  0.001 -0.002  0.021 -0.044  0.023
 [8,]  0.005 -0.002  0.020 -0.005  0.017  0.017  0.020  0.001 -0.021  0.019  0.023 -0.042
 [9,]  0.016  0.000  0.011  0.003  0.015  0.014 -0.041  0.008  0.033 -0.048  0.026  0.022
[10,]  0.020  0.000  0.008  0.001  0.009  0.010 -0.034  0.012  0.022 -0.059  0.030  0.029
[11,] -0.003  0.028  0.003  0.028  0.002  0.002 -0.019  0.014  0.005  0.028 -0.058  0.030
[12,]  0.041 -0.006  0.000  0.003  0.010  0.009 -0.007 -0.027  0.034 -0.052  0.027  0.025
[13,] -0.003  0.049 -0.007  0.028 -0.003 -0.002  0.015 -0.021  0.006  0.024 -0.057  0.033
[14,]  0.001 -0.006  0.043 -0.003  0.020  0.020  0.029 -0.021 -0.008  0.025  0.032 -0.057
[15,]  0.003  0.028 -0.004  0.049 -0.008 -0.009  0.008  0.010 -0.018  0.030 -0.061  0.030
[16,]  0.012  0.000  0.018 -0.008  0.035  0.010  0.028  0.015 -0.043  0.020  0.032 -0.052
[17,]  0.007  0.000  0.021 -0.007  0.014  0.040  0.021  0.010 -0.031  0.031  0.028 -0.060
[18,]  0.035  0.000  0.002  0.001  0.003  0.003 -0.008  0.001  0.007 -0.003 -0.001  0.004
[19,]  0.000  0.037  0.000  0.037  0.000  0.000  0.000  0.000  0.000 -0.001  0.003 -0.001
[20,]  0.002  0.000  0.035 -0.001  0.034  0.034  0.008 -0.001 -0.007  0.004 -0.002 -0.002
[21,]  0.035  0.000  0.002  0.001  0.003  0.003 -0.008  0.001  0.007 -0.003 -0.001  0.004
[22,]  0.000  0.037  0.000  0.037  0.000  0.000  0.000  0.000  0.000 -0.001  0.003 -0.001
[23,]  0.002  0.000  0.035 -0.001  0.034  0.034  0.008 -0.001 -0.007  0.004 -0.002 -0.002
[24,]  0.035  0.001  0.004  0.001  0.004  0.004 -0.013  0.005  0.008 -0.005 -0.001  0.006
[25,]  0.035  0.001  0.004  0.001  0.004  0.004 -0.013  0.004  0.009 -0.004 -0.001  0.006
[26,]  0.000  0.037  0.002  0.037  0.001  0.001 -0.006  0.004  0.002 -0.004  0.003  0.001
 [ getOption("max.print") 에 도달했습니다 -- 12 행들을 생략합니다 ]
> 
> diag(gi_LHS)
 [1] 0.71249489 0.63796043 0.48888695 0.48884202 0.48683355 0.48888695 0.48884202 0.48683355 0.43752765 0.42149993 0.43187307 0.43654180 0.43174939 0.42835526 0.42999170 0.43135400 0.42369129
[18] 0.07228834 0.07412064 0.07223240 0.07228834 0.07412064 0.07223240 0.06971005 0.06920679 0.07324710 0.07144132 0.07300684 0.07088560 0.07330810 0.06992155 0.06936149 0.20722490 0.19551685
[35] 0.21077837 0.21895820 0.21092827 0.21603916
> 
> # solution
> sol = gi_LHS %*% RHS
> 
> # fixed effect : sex effect
> sol_f1 = sol[1:2]
> sol_f1
[1] 6.004142 8.242687
> 
> # direct EBV
> sol_dbv = sol[3:17]
> 
> # associative EBV
> sol_sbv = sol [18:32]
> 
> sol_bv = cbind(sol_dbv, sol_sbv, sol_dbv + sol_sbv)
> sol_bv
         sol_dbv      sol_sbv            
 [1,]  0.2955862 -0.044464992  0.25112116
 [2,] -0.4833830  0.027828018 -0.45555496
 [3,]  0.1877968  0.016636974  0.20443380
 [4,]  0.2955862 -0.044464992  0.25112116
 [5,] -0.4833830  0.027828018 -0.45555496
 [6,]  0.1877968  0.016636974  0.20443380
 [7,]  0.1247823 -0.075976274  0.04880605
 [8,]  0.5218020 -0.098832417  0.42296958
 [9,] -0.8735821  0.008947170 -0.86463497
[10,]  0.5357603 -0.003051276  0.53270903
[11,] -0.4876411  0.083313528 -0.40432758
[12,]  0.3991916  0.059707501  0.45889906
[13,] -0.5723087  0.019051373 -0.55325729
[14,]  0.1530269  0.004742619  0.15776954
[15,]  0.1989688  0.002097776  0.20106659
> 
> # random effect : group effect
> sol_r1 = sol[33:35]
> sol_r1
[1] -0.26871657  0.35903353 -0.09031695
> 
> # random effect : common environmental effect
> sol_r2 = sol[36:38]
> sol_r2
[1]  0.3327777 -0.5153337  0.1825559
>

결과가 책과 아주 약간 다르다. blupf90으로 샐행한 결과와 더 가깝다.

저작자표시 (새창열림)

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

혈통 그루핑 (0)	2023.05.08
Single-step GBLUP(ssGBLUP) (0)	2021.01.05
GBLUP Model with residual polygenic effects (0)	2020.12.29
SNP-BLUP Model with Polygenic Effects(unweighted analysis) (0)	2020.12.24
Mixed Linear Model(GBLUP model) for Computing DGV (0)	2020.12.21

H 행렬, H 행렬의 역행렬 유도

투정이 2021. 1. 19. 20:55

2021. 1. 19. 20:55

H 행렬 유도

H 행렬의 역행렬 유도

Deriving H matrix

Deriving Inverse of H matrix

single-stepGBLUP의 MME에는 H 행렬의 역행렬이 등장한다. 마치 NRM과 GRM의 적당한 결합처럼 보인다. 아니 무식해서 그렇게 생각했다. 그러나 NRM과 GRM의 적당한 결합이 아니다.

다시 GBLUP을 생각해 보자. GBLUP은 SNP genotyping한 개체들의 GRM을 만들어서 개체들의 GEBV를 추정한다. 집단의 모든 개체들의 GEBV를 추정하려면 모든 개체를 SNP genotyping 하여야 한다. 그러나 한정된 예산 때문에 또는 이미 개체들이 죽어서 genotyping을 할 수 없다. 그러면 죽은 선조까지 포함하여 전체 집단에 대해서 GBLUP을 할 수는 없는 걸까? 이런 질문에서 출발하여 일부 유전체 분석을 한 개체들의 gene content를 이용하여 유전체 분석을 하지 않은 개체(non genotyped animals)들의 gene content를 추정하여 GBLUP을 하자는 것이 ssGBLUP이고 그때 H matrix를 이용하는 것이다.

유전체 분석을 한 개체는 유전체 분석 결과를 이용하고, 유전체 분석을 안 한 개체는 유전체 분석을 한 개체로부터 gene content를 추정하여 만든 GRM이 바로 H matrix이다.

genotyped animals가 주어졌을 때 non genotyped animals의 gene contents의 분포를 이용하는 것이므로 정규 분포의 조건부 분포에 대한 이해가 필요할 수도 있다. 키워드로 conditional normal distribution, mean vector of conditional normal distribution, covariance matrix of conditional normal distribution 등을 이용할 수 있다.

H matrix는 매우 복잡하나, H 역행렬은 매우 간단한 형태를 취하고 있다. 놀라운 자연의 법칙이고, 이걸 알아내는 과학자들도 대단하다. Simple is beautiful.

저작자표시 (새창열림)

'Animal Breeding > Genomic Selection' 카테고리의 다른 글

pregsf90으로 유전체 근교 계수 구하기 (0)	2023.09.14
서로 다른 Chip 버전을 이용한 친자감정 및 부모 찾기 (1)	2023.04.18
SNP Effect를 이용한 DGV(or GEBV) 추정 (0)	2020.03.05
PostGSF90을 이용한 GWAS (0)	2020.03.05
SNP Marker Effect 추정 (0)	2020.03.02

blupf90으로 single-step GBLUP Model 풀기

투정이 2021. 1. 5. 17:06

2021. 1. 5. 17:06

# Single-step GBLUP Model(ssGBLUP)

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

# Raphael Mrode

# Example 11.6 p192

간단한 설명은 아래 포스팅 참조

2021/01/05 - [Animal Breeding/R for Genetic Evaluation] - Single-step GBLUP(ssGBLUP)

Single-step GBLUP(ssGBLUP)

# Single-step GBLUP Model(ssGBLUP) # Linear Models for the Prediction of Animal Breeding Values, 3rd Edition. # Raphael Mrode # Example 11.6 p192 single step GBLUP(ssGBLUP)은 유전체 자..

bhpark.tistory.com

Data(ssgblup_data.txt)

13 0 0 1 558 9 0.001792115
14 0 0 1 722 13.4 0.001385042
15 13 4 1 300 12.7 0.003333333
16 15 2 1 73 15.4 0.01369863
17 15 5 1 52 5.9 0.019230769
18 14 6 1 87 7.7 0.011494253
19 14 9 1 64 10.2 0.015625
20 14 9 1 103 4.8 0.009708738
21 1 3 1 13 7.6 0.076923077
22 14 8 1 125 8.8 0.008

animal, sire, dam, general mean, edc, dyd fat, 1/edc

Pedigree(ssgblup_pedi.txt)

SNP genotype(ssgblup_snpgenotype.txt)

18 11010202210000000000000000000000000000000000000000
19 00110202200000000000000000000000000000000000000000
20 01100102200000000000000000000000000000000000000000
21 20000122120000000000000000000000000000000000000000
22 00011202000000000000000000000000000000000000000000
23 01100102210000000000000000000000000000000000000000
24 10001102000000000000000000000000000000000000000000
25 00011202100000000000000000000000000000000000000000
26 10110201000000000000000000000000000000000000000000

Renumf90 파라미터 파일(renumf90_ssgblup.par)

# Parameter file for program renf90; it is translated to parameter
# file for BLUPF90 family programs.
DATAFILE
ssgblup_data.txt
TRAITS
6
FIELDS_PASSED TO OUTPUT
 
WEIGHT(S)
 7
RESIDUAL_VARIANCE
245
EFFECT
4 cross alpha
EFFECT
1 cross alpha
RANDOM
animal
FILE
ssgblup_pedi.txt
FILE_POS
1 2 3
SNP_FILE
ssgblup_snpgenotype.txt
PED_DEPTH
0
(CO)VARIANCES
35.241
OPTION no_quality_control
OPTION AlphaBeta 0.95 0.05
OPTION tunedG 0
OPTION thrStopCorAG 0.10
OPTION solv_method FSPAK

Renumf90 실행

renumf90 renumf90_ssgblup.par | tee renumf90_ssgblup.log

Renumf90 실행 화면

Renumf90 로그

 RENUMF90 version 1.145
 renumf90_ssgblup.par
 datafile:ssgblup_data.txt
 traits:           6
 R
   245.0    

 Processing effect  1 of type cross     
 item_kind=alpha     

 Processing effect  2 of type cross     
 item_kind=alpha     
 pedigree file name  "ssgblup_pedi.txt"
 positions of animal, sire, dam, alternate dam, yob, and group     1     2     3     0     0     0     0
 SNP file name  "ssgblup_snpgenotype.txt"
 all pedigrees to be included
 Reading (CO)VARIANCES:           1 x           1

 Maximum size of character fields: 20

 Maximum size of record (max_string_readline): 800

 Maximum number of fields for input file (max_field_readline): 100

 Pedigree search method (ped_search): convention

 Order of pedigree animals (animal_order): default

 Order of UPG (upg_order): default

 Missing observation code (missing): 0

 hash tables for effects set up
 first 3 lines of the data file (up to 70 characters)
    13 0 0 1 558 9 0.001792115
    14 0 0 1 722 13.4 0.001385042
    15 13 4 1 300 12.7 0.003333333
 read           10  records
 table with            1  elements sorted
 added count
 Effect group            1  of column            1  with            1  levels
 table expanded from        10000  to        10000  records
 added count
 Effect group            2  of column            1  with           10  levels
 wrote statistics in file "renf90.tables"

 Basic statistics for input data  (missing value code is '0')
 Pos  Min         Max         Mean        SD                 N
   6    4.8000      15.400      9.5500      3.3890          10

 random effect with SNPs  2
 type: animal    
 file: ssgblup_snpgenotype.txt
 read SNPs           9  records
 Effect group            2  of column            1  with           14  levels

 random effect   2
 type:animal    
 opened output pedigree file "renadd02.ped"
 read           26  pedigree records
 loaded           12  parent(s) in round            0

 Pedigree checks
 
 Number of animals with records                  =           10
 Number of animals with genotypes                =            9
 Number of animals with records or genotypes     =           14
 Number of animals with genotypes and no records =            4
 Number of parents without records or genotypes  =           12
 Total number of animals                         =           26

 Wrote cross reference IDs for SNP file "ssgblup_snpgenotype.txt_XrefID"

 Wrote parameter file "renf90.par"
 Wrote renumbered data "renf90.dat" 10 records

Renumf90 실행 결과 생성된 파일

renadd02.ped

14 4 26 1 0 12 0 0 0 26
1 0 0 3 0 0 1 1 0 13
2 15 17 1 0 12 1 0 0 21
19 0 0 3 0 0 0 0 1 5
3 4 23 1 0 12 1 0 0 19
4 0 0 3 0 0 1 8 0 14
5 4 22 1 0 12 1 0 0 22
17 0 0 3 0 0 0 0 1 3
22 0 0 3 0 0 0 0 1 8
6 1 18 1 0 2 1 2 0 15
11 4 25 1 0 12 0 0 0 23
24 0 0 3 0 0 0 0 1 10
15 0 0 3 0 0 0 1 0 1
20 0 0 3 0 0 0 0 1 6
7 6 16 1 0 2 1 0 0 16
12 4 24 1 0 12 0 0 0 24
25 0 0 3 0 0 0 0 1 11
18 0 0 3 0 0 0 0 1 4
8 6 19 1 0 2 1 0 0 17
23 0 0 3 0 0 0 0 2 9
13 4 21 1 0 12 0 0 0 25
26 0 0 3 0 0 0 0 1 12
9 4 23 1 0 12 1 0 0 20
16 0 0 3 0 0 0 0 1 2
10 4 20 1 0 12 1 0 0 18
21 0 0 3 0 0 0 0 1 7

Cross Reference ID(ssgblup_snpgenotype.txt_XrefID) - genotype을 가지고 있는 개체들의 renumbered ID

renf90.dat

 9 0.001792115 1 1
 13.4 0.001385042 1 4
 12.7 0.003333333 1 6
 15.4 0.01369863 1 7
 5.9 0.019230769 1 8
 7.7 0.011494253 1 10
 10.2 0.015625 1 3
 4.8 0.009708738 1 9
 7.6 0.076923077 1 2
 8.8 0.008 1 5

renf90.par

# BLUPF90 parameter file created by RENUMF90
DATAFILE
 renf90.dat
NUMBER_OF_TRAITS
           1
NUMBER_OF_EFFECTS
           2
OBSERVATION(S)
    1
WEIGHT(S)
           2
EFFECTS: POSITIONS_IN_DATAFILE NUMBER_OF_LEVELS TYPE_OF_EFFECT[EFFECT NESTED]
 3         1 cross 
 4        26 cross 
RANDOM_RESIDUAL VALUES
   245.00    
 RANDOM_GROUP
     2
 RANDOM_TYPE
 add_animal   
 FILE
renadd02.ped                                                                    
(CO)VARIANCES
   35.241    
OPTION SNP_file ssgblup_snpgenotype.txt
OPTION no_quality_control
OPTION AlphaBeta 0.95 0.05
OPTION tunedG 0
OPTION thrStopCorAG 0.10
OPTION solv_method FSPAK

본 예제의 genotype은 제대로 된 genotype이 아니다. 그래서 다음 프로그램이 정상적으로 동작하지 않는다. 그래서 no_quality_control, thrStopCorAG 등의 옵션을 넣어서 강제로 실행을 하는 것이다. 실제 분석에서는 세심한 주의를 기울여서 옵션을 주고 문제점이 있나 없나 확인을 해야 한다. 그렇지 않으면 유전체 자료를 이용하는 것이 독이 될 수도 있다.

중간 과정을 점검하기 위하여 H 역행렬, A 역행렬, GimA22i 등등이 필요하면 위의 "OPTION solv_method FSPAK"을 지우고 다음을 추가하여 pregsf90을 실행한다.

OPTION saveAscii
OPTION saveHinv
OPTION saveAinv
OPTION saveHinvOrig
OPTION saveAinvOrig
OPTION saveDiagGOrig
OPTION saveGOrig
OPTION saveA22Orig
OPTION saveGimA22iOrig

여기서는 실행하지 않는다.

실제에서는 pregsf90을 실행하여 유전체 자료의 문제점이 있는지 없는지 반드시 확인하여야 한다. pregsf90은 문제점을 찾아내고 수정하기 위한 많은 옵션을 제공하는데 반드시 의미하는 바가 무엇인지 확인하고 진행하기를 바란다. 그리고 pregsf90, blupf90은 문제가 있어도 대충 덮고 실행을 계속하게 한다. 그러나 그렇게 하지 말기를 바란다. 문제가 발견되면 문제를 일으킨 유전체 자료가 무엇인지 확인하고 분석비용이 아깝긴 하지만 넣는 것이 바람직한 것인지 아니면 버리는 것이 나은 것인지 판단을 하고 진행해야 한다. 보통 버리는 것이 정신 건강상 좋다. 유전체 자료를 넣는다고 무조건 정확도가 올라가서 개량량이 증가하는 것은 아니다. 양질의 유전체 자료를 넣을 때만 그렇게 된다. 당연한 얘기가 실제에선 무시되고 마구 분석하는 현실이 안타깝기만 하다. 모든 것이 그렇듯이 유전체 자료를 이용한 유전체 선발도 마법의 도구가 아니다.

blupf90 실행하기 - renumf90이 만들어낸 파라미터 파일을 이용하여 blupf90 실행

다음 명령어로 실행

blupf90 renf90.par | tee blupf90_ssgblup.log

실행 화면

실행 로그

renf90.par
     BLUPF90 ver. 1.68

 Parameter file:             renf90.par
 Data file:                  renf90.dat
 Number of Traits             1
 Number of Effects            2
 Position of Observations      1
 Position of Weight (1)        2
 Value of Missing Trait/Observation           0

EFFECTS
 #  type                position (2)        levels   [positions for nested]
    1  cross-classified       3         1
    2  cross-classified       4        26

 Residual (co)variance Matrix
  245.00    

 Random Effect(s)    2
 Type of Random Effect:      additive animal
 Pedigree File:              renadd02.ped                                                                                                                                                                                                                                              
 trait   effect    (CO)VARIANCES
  1       2     35.24    

 REMARKS
  (1) Weight position 0 means no weights utilized
  (2) Effect positions of 0 for some effects and traits means that such
      effects are missing for specified traits
 

 * The limited number of OpenMP threads = 4

 * solving method (default=PCG):FSPAK               
 

 Options read from parameter file for genomic
 * SNP format: BLUPF90 standard (text)
 * SNP file: ssgblup_snpgenotype.txt
 * SNP Xref file: ssgblup_snpgenotype.txt_XrefID
 * No Quality Control Checks !!!!! (default .false.):  T
 * Set the threshold to Stop the analysis if low cor(A22,G) (default 0.3):  0.1
 * Create a tuned G (default = 2): 0
 * AlphaBeta defaults alpha=0.95, beta=0.05) :  0.95  0.05
 Data record length =            4
 # equations =           27
 G
  35.241    
 read           10  records in   0.9531250      s,                      21 
  nonzeroes
  read           26  additive pedigrees
 
 *--------------------------------------------------------------*
 *                 Genomic Library: Dist Version 1.281          *
 *                                                              *
 *             Optimized OpenMP Version -  4 threads            *
 *                                                              *
 *  Modified relationship matrix (H) created for effect:   2    *
 *--------------------------------------------------------------*
 
 Read 26 animals from pedigree file: "D:\users\bhpark\2021\job\공부하기\01_Linear Models for the Prediction of Animal Breeding Values_3rd_Mrode\26_ssGBLUP_blupf90\renadd02.ped"
 Number of Genotyped Animals: 9

 Creating A22 
    Extracting subset of: 19 pedigrees from: 26 elapsed time:     0.0000
    Calculating A22 Matrix by Colleau OpenMP...elapsed time: .0000
    Numbers of threads=1 4

 Reading SNP file
    Column position in file for the first marker: 4
    Format to read SNP file: (3x,400000i1)                                     
    Number of SNPs: 50
    Format: integer genotypes (0 to 5) to double-precision array
    Number of Genotyped animals: 9
    Reading SNP file elapsed time: .00
 
 Statistics of alleles frequencies in the current population
    N:             50
    Mean:       0.074
    Min:        0.000
    Max:        0.944
    Var:        0.038
 
 
 Quality Control - Monomorphic SNPs Exist - NOT REMOVED: 40

 Genotypes missings (%):  0.000

 Calculating G Matrix 
    Dgemm MKL #threads=     1    4 Elapsed omp_get_time:     0.0000
 
 Scale by Sum(2pq). Average:   3.19135802469136     

 Detecting samples with similar genotypes
    elapsed time=     0.0
 
 Blend G as alpha*G + beta*A22: (alpha,beta)     0.950     0.050

 Frequency - Diagonal of G
    N:           9
    Mean:        1.057
    Min:         0.565
    Max:         2.648
    Range:       0.104
    Class:     20
 
  #Class       Class   Count
       1  0.5645           1
       2  0.6687           1
       3  0.7729           1
       4  0.8771           4
       5  0.9813           1
       6   1.085           0
       7   1.190           0
       8   1.294           0
       9   1.398           0
      10   1.502           0
      11   1.606           0
      12   1.711           0
      13   1.815           0
      14   1.919           0
      15   2.023           0
      16   2.127           0
      17   2.232           0
      18   2.336           0
      19   2.440           0
      20   2.544           1
      21   2.648           0
 

 Check for diagonal of genomic relationship matrix

 ** High Diagonal of genotype 4  2.65 Not Removed
 ** Low Diagonal of genotype 8  0.56 Not Removed

 Check for diagonal of genomic relationship matrix, genotypes not removed: 2

 ------------------------------
  Final Pedigree-Based Matrix 
 ------------------------------
 
 Statistic of Rel. Matrix A22
                            N      Mean       Min       Max       Var
     Diagonal               9     1.000     1.000     1.000     0.000
     Off-diagonal          72     0.201     0.000     0.500     0.013
 

 ----------------------
  Final Genomic Matrix 
 ----------------------
 
 Statistic of Genomic Matrix
                            N      Mean       Min       Max       Var
     Diagonal               9     1.057     0.565     2.648     0.377
     Off-diagonal          72    -0.116    -0.742     0.725     0.142
 

 Correlation of Genomic Inbreeding and Pedigree Inbreeding
     Variance of Y is Zero !!
     Correlation:     0.0000
 
 Diagonal elements
    Estimating Regression Coefficients G = b0 11' + b1 A + e
    Regression coefficients b0 b1 =        NaN       NaN

    Correlation diagonal elements G & A       NaN
 
 All elements - Diagonal / Off-Diagonal
    Estimating Regression Coefficients G = b0 11' + b1 A + e
    Regression coefficients b0 b1 =     -0.395     1.413

    Correlation all elements G & A     0.708
 
 Off-Diagonal
    Using 56 elements from A22 >= .02000

    Estimating Regression Coefficients G = b0 11' + b1 A + e
    Regression coefficients b0 b1 =     -0.483     1.647

    Correlation Off-Diagonal elements G & A     0.212
 
 ***********************************************************************
 * CORRELATION FOR OFF-DIAGONALS G & A22 IS LOW THAN  0.50  !!!!!  *
 * MISIDENTIFIED GENOMIC SAMPLES OR POOR QUALITY GENOMIC DATA *
 ***********************************************************************

 Creating A22-inverse 
    Inverse LAPACK MKL dpotrf/i  #threads=    1    4 Elapsed omp_get_time:     0.0000

 ----------------------
  Final A22 Inv Matrix 
 ----------------------
 
 Statistic of Inv. Rel. Matrix A22
                            N      Mean       Min       Max       Var
     Diagonal               9     1.233     1.000     1.429     0.017
     Off-diagonal          72    -0.101    -0.571     0.000     0.009
 
 
 Creating G-inverse 
    Inverse LAPACK MKL dpotrf/i  #threads=    1    4 Elapsed omp_get_time:     0.0000

 --------------------------
  Final Genomic Inv Matrix 
 --------------------------
 
 Statistic of Inv. Genomic Matrix
                            N      Mean       Min       Max       Var
     Diagonal               9     5.227     1.593    11.227     9.161
     Off-diagonal          72     0.308    -5.482     3.166     3.142
 

 Check for diagonal of Inverse Genomic - Inverse of pedigree relationship matrix


 ------------------------------
  Final G Inv - A22 Inv Matrix 
 ------------------------------
 
 Statistic of Inv. Genomic - A22 Matrix
                            N      Mean       Min       Max       Var
     Diagonal               9     3.994     0.593    10.021     8.878
     Off-diagonal          72     0.408    -5.355     3.293     3.069
 

*--------------------------------------------------*
* Setup Genomic Done !!!, elapsed time:      0.049 *
*--------------------------------------------------*

 finished peds in    1.156250      s,                     108  nonzeroes
 solutions stored in file: "solutions"

blupf90 실행결과 생기는 파일

sum2pq

   3.19135802469136

solutions

trait/effect level  solution
   1   1         1          8.39032676
   1   2         1          0.00291425
   1   2         2         -0.02384446
   1   2         3          0.00746559
   1   2         4          0.00433541
   1   2         5          0.00839707
   1   2         6          0.00559388
   1   2         7          0.01313654
   1   2         8         -0.00236465
   1   2         9          0.00100295
   1   2        10         -0.00244449
   1   2        11         -0.00298365
   1   2        12         -0.00082301
   1   2        13          0.00793245
   1   2        14          0.00483534
   1   2        15         -0.01192223
   1   2        16          0.00689306
   1   2        17         -0.01192223
   1   2        18          0.00275784
   1   2        19         -0.00344106
   1   2        20         -0.00307480
   1   2        21          0.00384316
   1   2        22          0.00415291
   1   2        23          0.00206656
   1   2        24         -0.00199381
   1   2        25         -0.00343423
   1   2        26          0.00177842

original ID 로 정렬한 개체의 gebv

rid	oid	sol
15	1	-0.01192
16	2	0.00689
17	3	-0.01192
18	4	0.00276
19	5	-0.00344
20	6	-0.00307
21	7	0.00384
22	8	0.00415
23	9	0.00207
24	10	-0.00199
25	11	-0.00343
26	12	0.00178
1	13	0.00291
4	14	0.00434
6	15	0.00559
7	16	0.01314
8	17	-0.00236
10	18	-0.00244
3	19	0.00747
9	20	0.00100
2	21	-0.02384
5	22	0.00840
11	23	-0.00298
12	24	-0.00082
13	25	0.00793
14	26	0.00484

해가 책과도 다르고 blupf90 팀의 Yutaka Masuda의 해와도 다르다. 위의 참고 포스팅하고는 같다. 이유는 참고 포스팅을 보기 바란다.

저작자표시 (새창열림)

'Animal Breeding > BLUPF90' 카테고리의 다른 글

inbupgf90을 이용한 근교계수 계산과 교배계획 세우기 (0)	2023.05.25
blupf90으로 social interaction model 풀기 (1)	2021.04.21
blupf90으로 GBLUP Model with residual polygenic effects 풀기 (0)	2020.12.29
blupf90으로 gblup model 풀기 (0)	2020.12.29
blupf90으로 SNP-BLUP Model(with polygenic effects unweighted analysis) 풀기 (0)	2020.12.24

Single-step GBLUP(ssGBLUP)

투정이 2021. 1. 5. 15:57

2021. 1. 5. 15:57

# Single-step GBLUP Model(ssGBLUP)

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

# Raphael Mrode

# Example 11.6 p192

single step GBLUP(ssGBLUP)은 유전체 자료를 가지고 있는 개체와 가지고 있지 않은 개체를 동시에 분석할 수 있는 방법이다. 또한 SNP efffect를 구한 다음 DGV 또는 GEBV를 구하는 것이 아니라 SNP effect를 구하는 단계를 생략하고 한 번에 유전체 자료를 가지고 있는 개체는 GEBV, 없는 개체는 EBV를 구한다.

ssGBLUP에 포함하는 개체는 다음의 세 가지로 분류할 수 있다.

1. 능력검정 기록만 가지고 있는 개체 - 유전체 자료를 생성하기 전에, 즉 과거에 유전평가에 포함된 개체

2. 능력검정 기록과 유전체 자료를 가지고 있는 개체 - 소위 참조집단(reference population)이라고 하는 개체들로 GEBV를 추정한다. 참조집단이 많아야 유전체 선발(genomic selection)의 효율(아래 3번 개체들의 정확도 향상)이 높아진다.

3. 유전체 자료만 가지고 있는 개체 - 보통 어린 개체들로 아직 능력검정을 할 단계에 이르지 못한 개체. 이들에 대해 GEBV를 계산하고 조기 선발함으로써 세대간격(generation interval)을 줄여 유전적 개량량을 높일 수 있다. 이들은 시간이 지나 능력검정할 단계에 다다르면 능력검정을 하여 표현형 자료를 확보하여 참조집단이 된다.

모형 또는 MME는 위의 책 또는 아래 포스팅을 참고한다.

masuday.github.io/blupf90_tutorial/mrode_c11ex116_ssgblup.html

Numerical examples from Mrode (2014)

Numerical examples from Mrode (2014) Yutaka Masuda September 2019 Back to index.html. Single-step approach Model We consider the standard animal model \[ \mathbf{y}=\mathbf{Xb}+\mathbf{Wu}+\mathbf{e}. \] If some animals are genotyped, their additive relati

masuday.github.io

Data(ssgblup_data.txt)

animal sire dam mean edc fat_dyd snp1 snp2 snp3 snp4 snp5 snp6 snp7 snp8 snp9 snp10
13 0 0 1 558 9 2 0 1 1 0 0 0 2 1 2
14 0 0 1 722 13.4 1 0 0 0 0 2 0 2 1 0
15 13 4 1 300 12.7 1 1 2 1 1 0 0 2 1 2
16 15 2 1 73 15.4 0 0 2 1 0 1 0 2 2 1
17 15 5 1 52 5.9 0 1 1 2 0 0 0 2 1 2
18 14 6 1 87 7.7 1 1 0 1 0 2 0 2 2 1
19 14 9 1 64 10.2 0 0 1 1 0 2 0 2 2 0
20 14 9 1 103 4.8 0 1 1 0 0 1 0 2 2 0
21 1 3 1 13 7.6 2 0 0 0 0 1 2 2 1 2
22 14 8 1 125 8.8 0 0 0 1 1 2 0 2 0 0
23 14 11 1 93 9.8 0 1 1 0 0 1 0 2 2 1
24 14 10 1 66 9.2 1 0 0 0 1 1 0 2 0 0
25 14 7 1 75 11.5 0 0 0 1 1 2 0 2 1 0
26 14 12 1 33 13.3 1 0 1 1 0 2 0 1 0 0

Pedigree(ssgblup_pedi.txt)

animal sire dam
1 0 0
2 0 0
3 0 0
4 0 0
5 0 0
6 0 0
7 0 0
8 0 0
9 0 0
10 0 0
11 0 0
12 0 0
13 0 0
14 0 0
15 13 4
16 15 2
17 15 5
18 14 6
19 14 9
20 14 9
21 1 3
22 14 8
23 14 11
24 14 10
25 14 7
26 14 12

R 프로그램

# single step GBLUP Model

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

# Raphael Mrode

# Example 11.6 p192

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

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

# function to make numerator relationship matrix
makenrm = function(pedi) {
  n = nrow(pedi)
  nrm = matrix(c(0), nrow = n, ncol = n)
  
  for (i in 1:n) {
    animal = pedi[i, 1]
    sire = pedi[i, 2]
    dam = pedi[i, 3]
    
    if (sire == 0 & dam == 0) {
      # both parents unknown
      nrm[animal, animal] = 1
    } else if (sire != 0 & dam == 0) {
      # sire known
      for (j in 1:animal - 1) {
        nrm[j, animal] = 0.5 * (nrm[j, sire])
        nrm[animal, j] = nrm[j, animal]
      }
      nrm[animal, animal] = 1
    } else if (sire == 0 & dam != 0) {
      # dam known
      for (j in 1:animal - 1) {
        nrm[j, animal] = 0.5 * (nrm[j, dam])
        nrm[animal, j] = nrm[j, animal]
      }
      nrm[animal, animal] = 1
    } else {
      # both parents known
      for (j in 1:animal - 1) {
        nrm[j, animal] = 0.5 * (nrm[j, sire] + nrm[j, dam])
        nrm[animal, j] = nrm[j, animal]
      }
      nrm[animal, animal] = 1 + 0.5 * nrm[sire, dam]
    }
  }
  return(nrm)
}

# set working directory 
setwd(".") 

# print working directory 
getwd()

# read data
data = read.table("ssgblup_data.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
data

# read pedigree : animal sire dam
pedi = read.table("ssgblup_pedi.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
pedi

# variance component and ratio
sigma_a = 35.241
sigma_e = 245
sigma_a
sigma_e

# variance ratio, alpha
alpha = sigma_e / sigma_a
alpha

# SNP 
M_all = as.matrix(data[6:14, 7:16])
M_all

#  mean of each SNP(allele frequency)
snp_mean = colMeans(M_all) / 2
snp_mean

#  sum of 2pq(heterozygote ratio)
sum2pq = sum(2 * snp_mean * (1 - snp_mean))
sum2pq

# observation
y = data[1:10, 6]
y

# design matrix of fixed effect(general mean)
X = matrix(rep(1, 10), 10, 1)
X

# design matrix for animal(parents : 12, animals with record : 10, animals without record : 4)
W = cbind(matrix(c(0), 10, 12), diag(10), matrix(c(0), 10, 4))
W

# genomic relationship matrix, G
Z = sweep(M_all, 2, snp_mean * 2)
round(Z, 3)
G = Z %*% t(Z) / sum2pq
round(G, 3)

# Numerator Relationship Matrix
A = makenrm(pedi)
A
A22 = A[18:26, 18:26]
A22

# weighted G
wt = 0.95
Gw = wt * G + (1 - wt) * A22
round(Gw, 3)
Gwi = ginv(Gw)
round(Gwi, 3)

# inverse matrix of NRM 
Ai = ainv(pedi) 
Ai

# H matrix
H0 = matrix(c(0), 26, 26)
H0[18:26, 18:26] = Gwi - ginv(A22)
H0
Hi = Ai + H0
round(Hi, 3)

# weight R
R = diag(data[1:10, 5])
R
Ri = ginv(R)
Ri

# Left Hand Side
LHS = rbind(
  cbind(t(X) %*% Ri %*% X, t(X) %*% Ri %*% W), 
  cbind(t(W) %*% Ri %*% X, t(W) %*% Ri %*% W + Hi * alpha))
round(LHS, 3)

# Right Hand Side
RHS = rbind(t(X) %*% Ri %*% y, 
            t(W) %*% Ri %*% y)
round(RHS, 3)

# Solutions

# generalized inverse of LHS
gi_LHS = ginv(LHS)
round(gi_LHS, 3)

# solution
sol = gi_LHS %*% RHS
round(sol, 6)

# general mean effect
gme = sol[1]
gme

# (g)ebv of animals
ebv = sol[2:27]
ebv

RStudio에서 실행하는 화면

실행 결과

> # single step GBLUP Model
> 
> # Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.
> 
> # Raphael Mrode
> 
> # Example 11.6 p192
> 
> # MASS packages
> if (!("MASS" %in% installed.packages())) {
+   install.packages("MASS", dependencies = TRUE)
+ }
> library(MASS)
> 
> # function to make inverse of numerator relationship matrix
> ainv = function(pedi) {
+   n = nrow(pedi)
+   Ainv = matrix(c(0), nrow = n, ncol = n)
+   
+   for (i in 1:n) {
+     animal = pedi[i, 1]
+     sire = pedi[i, 2]
+     dam = pedi[i, 3]
+     
+     if (sire == 0 & dam == 0) {
+       # both parents unknown
+       alpha = 1
+       Ainv[animal, animal] = alpha + Ainv[animal, animal]
+     } else if (sire != 0 & dam == 0) {
+       # sire known
+       alpha = 4/3
+       Ainv[animal, animal] = alpha + Ainv[animal, animal]
+       Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
+       Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
+       Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
+     } else if (sire == 0 & dam != 0) {
+       # dam known
+       alpha = 4/3
+       Ainv[animal, animal] = alpha + Ainv[animal, animal]
+       Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
+       Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
+       Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
+     } else {
+       # both parents known
+       alpha = 2
+       Ainv[animal, animal] = alpha + Ainv[animal, animal]
+       Ainv[animal, sire] = -alpha/2 + Ainv[animal, sire]
+       Ainv[sire, animal] = -alpha/2 + Ainv[sire, animal]
+       Ainv[animal, dam] = -alpha/2 + Ainv[animal, dam]
+       Ainv[dam, animal] = -alpha/2 + Ainv[dam, animal]
+       Ainv[sire, sire] = alpha/4 + Ainv[sire, sire]
+       Ainv[sire, dam] = alpha/4 + Ainv[sire, dam]
+       Ainv[dam, sire] = alpha/4 + Ainv[dam, sire]
+       Ainv[dam, dam] = alpha/4 + Ainv[dam, dam]
+     }
+   }
+   return(Ainv)
+ }
> 
> # function to make numerator relationship matrix
> makenrm = function(pedi) {
+   n = nrow(pedi)
+   nrm = matrix(c(0), nrow = n, ncol = n)
+   
+   for (i in 1:n) {
+     animal = pedi[i, 1]
+     sire = pedi[i, 2]
+     dam = pedi[i, 3]
+     
+     if (sire == 0 & dam == 0) {
+       # both parents unknown
+       nrm[animal, animal] = 1
+     } else if (sire != 0 & dam == 0) {
+       # sire known
+       for (j in 1:animal - 1) {
+         nrm[j, animal] = 0.5 * (nrm[j, sire])
+         nrm[animal, j] = nrm[j, animal]
+       }
+       nrm[animal, animal] = 1
+     } else if (sire == 0 & dam != 0) {
+       # dam known
+       for (j in 1:animal - 1) {
+         nrm[j, animal] = 0.5 * (nrm[j, dam])
+         nrm[animal, j] = nrm[j, animal]
+       }
+       nrm[animal, animal] = 1
+     } else {
+       # both parents known
+       for (j in 1:animal - 1) {
+         nrm[j, animal] = 0.5 * (nrm[j, sire] + nrm[j, dam])
+         nrm[animal, j] = nrm[j, animal]
+       }
+       nrm[animal, animal] = 1 + 0.5 * nrm[sire, dam]
+     }
+   }
+   return(nrm)
+ }
> 
> # set working directory 
> setwd(".") 
> 
> # print working directory 
> getwd()
[1] "D:/users/bhpark/2021/job/공부하기/01_Linear Models for the Prediction of Animal Breeding Values_3rd_Mrode/25_ssGBLUP"
> 
> # read data
> data = read.table("ssgblup_data.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
> data
   animal sire dam mean edc fat_dyd snp1 snp2 snp3 snp4 snp5 snp6 snp7 snp8 snp9 snp10
1      13    0   0    1 558     9.0    2    0    1    1    0    0    0    2    1     2
2      14    0   0    1 722    13.4    1    0    0    0    0    2    0    2    1     0
3      15   13   4    1 300    12.7    1    1    2    1    1    0    0    2    1     2
4      16   15   2    1  73    15.4    0    0    2    1    0    1    0    2    2     1
5      17   15   5    1  52     5.9    0    1    1    2    0    0    0    2    1     2
6      18   14   6    1  87     7.7    1    1    0    1    0    2    0    2    2     1
7      19   14   9    1  64    10.2    0    0    1    1    0    2    0    2    2     0
8      20   14   9    1 103     4.8    0    1    1    0    0    1    0    2    2     0
9      21    1   3    1  13     7.6    2    0    0    0    0    1    2    2    1     2
10     22   14   8    1 125     8.8    0    0    0    1    1    2    0    2    0     0
11     23   14  11    1  93     9.8    0    1    1    0    0    1    0    2    2     1
12     24   14  10    1  66     9.2    1    0    0    0    1    1    0    2    0     0
13     25   14   7    1  75    11.5    0    0    0    1    1    2    0    2    1     0
14     26   14  12    1  33    13.3    1    0    1    1    0    2    0    1    0     0
> 
> # read pedigree : animal sire dam
> pedi = read.table("ssgblup_pedi.txt", header = TRUE, sep = " ", stringsAsFactors = FALSE)
> pedi
   animal sire dam
1       1    0   0
2       2    0   0
3       3    0   0
4       4    0   0
5       5    0   0
6       6    0   0
7       7    0   0
8       8    0   0
9       9    0   0
10     10    0   0
11     11    0   0
12     12    0   0
13     13    0   0
14     14    0   0
15     15   13   4
16     16   15   2
17     17   15   5
18     18   14   6
19     19   14   9
20     20   14   9
21     21    1   3
22     22   14   8
23     23   14  11
24     24   14  10
25     25   14   7
26     26   14  12
> 
> # variance component and ratio
> sigma_a = 35.241
> sigma_e = 245
> sigma_a
[1] 35.241
> sigma_e
[1] 245
> 
> # variance ratio, alpha
> alpha = sigma_e / sigma_a
> alpha
[1] 6.95213
> 
> # SNP 
> M_all = as.matrix(data[6:14, 7:16])
> M_all
   snp1 snp2 snp3 snp4 snp5 snp6 snp7 snp8 snp9 snp10
6     1    1    0    1    0    2    0    2    2     1
7     0    0    1    1    0    2    0    2    2     0
8     0    1    1    0    0    1    0    2    2     0
9     2    0    0    0    0    1    2    2    1     2
10    0    0    0    1    1    2    0    2    0     0
11    0    1    1    0    0    1    0    2    2     1
12    1    0    0    0    1    1    0    2    0     0
13    0    0    0    1    1    2    0    2    1     0
14    1    0    1    1    0    2    0    1    0     0
> 
> #  mean of each SNP(allele frequency)
> snp_mean = colMeans(M_all) / 2
> snp_mean
     snp1      snp2      snp3      snp4      snp5      snp6      snp7      snp8      snp9     snp10 
0.2777778 0.1666667 0.2222222 0.2777778 0.1666667 0.7777778 0.1111111 0.9444444 0.5555556 0.2222222 
> 
> #  sum of 2pq(heterozygote ratio)
> sum2pq = sum(2 * snp_mean * (1 - snp_mean))
> sum2pq
[1] 3.191358
> 
> # observation
> y = data[1:10, 6]
> y
 [1]  9.0 13.4 12.7 15.4  5.9  7.7 10.2  4.8  7.6  8.8
> 
> # design matrix of fixed effect(general mean)
> X = matrix(rep(1, 10), 10, 1)
> X
      [,1]
 [1,]    1
 [2,]    1
 [3,]    1
 [4,]    1
 [5,]    1
 [6,]    1
 [7,]    1
 [8,]    1
 [9,]    1
[10,]    1
> 
> # design matrix for animal(parents : 12, animals with record : 10, animals without record : 4)
> W = cbind(matrix(c(0), 10, 12), diag(10), matrix(c(0), 10, 4))
> W
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
 [1,]    0    0    0    0    0    0    0    0    0     0     0     0     1     0     0     0     0     0     0     0     0     0     0     0     0     0
 [2,]    0    0    0    0    0    0    0    0    0     0     0     0     0     1     0     0     0     0     0     0     0     0     0     0     0     0
 [3,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     1     0     0     0     0     0     0     0     0     0     0     0
 [4,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     1     0     0     0     0     0     0     0     0     0     0
 [5,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     1     0     0     0     0     0     0     0     0     0
 [6,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     1     0     0     0     0     0     0     0     0
 [7,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     1     0     0     0     0     0     0     0
 [8,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     1     0     0     0     0     0     0
 [9,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0     1     0     0     0     0     0
[10,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0     0     0     0     0     1     0     0     0     0
> 
> # genomic relationship matrix, G
> Z = sweep(M_all, 2, snp_mean * 2)
> round(Z, 3)
     snp1   snp2   snp3   snp4   snp5   snp6   snp7   snp8   snp9  snp10
6   0.444  0.667 -0.444  0.444 -0.333  0.444 -0.222  0.111  0.889  0.556
7  -0.556 -0.333  0.556  0.444 -0.333  0.444 -0.222  0.111  0.889 -0.444
8  -0.556  0.667  0.556 -0.556 -0.333 -0.556 -0.222  0.111  0.889 -0.444
9   1.444 -0.333 -0.444 -0.556 -0.333 -0.556  1.778  0.111 -0.111  1.556
10 -0.556 -0.333 -0.444  0.444  0.667  0.444 -0.222  0.111 -1.111 -0.444
11 -0.556  0.667  0.556 -0.556 -0.333 -0.556 -0.222  0.111  0.889  0.556
12  0.444 -0.333 -0.444 -0.556  0.667 -0.556 -0.222  0.111 -1.111 -0.444
13 -0.556 -0.333 -0.444  0.444  0.667  0.444 -0.222  0.111 -0.111 -0.444
14  0.444 -0.333  0.556  0.444 -0.333  0.444 -0.222 -0.889 -1.111 -0.444
> G = Z %*% t(Z) / sum2pq
> round(G, 3)
        6      7      8      9     10     11     12     13     14
6   0.785  0.124  0.054  0.193 -0.398  0.228 -0.538 -0.120 -0.329
7   0.124  0.716  0.333 -0.781 -0.120  0.193 -0.573  0.159 -0.050
8   0.054  0.333  0.890 -0.538 -0.503  0.750 -0.329 -0.224 -0.433
9   0.193 -0.781 -0.538  2.735 -0.677 -0.050  0.124 -0.712 -0.294
10 -0.398 -0.120 -0.503 -0.677  0.925 -0.642  0.472  0.576  0.368
11  0.228  0.193  0.750 -0.050 -0.642  0.925 -0.468 -0.364 -0.573
12 -0.538 -0.573 -0.329  0.124  0.472 -0.468  0.959  0.124  0.228
13 -0.120  0.159 -0.224 -0.712  0.576 -0.364  0.124  0.542  0.019
14 -0.329 -0.050 -0.433 -0.294  0.368 -0.573  0.228  0.019  1.064
> 
> # Numerator Relationship Matrix
> A = makenrm(pedi)
> A
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
 [1,]  1.0  0.0  0.0 0.00  0.0  0.0  0.0  0.0  0.0   0.0   0.0   0.0  0.00   0.0   0.0  0.00  0.00  0.00  0.00  0.00   0.5  0.00  0.00  0.00  0.00  0.00
 [2,]  0.0  1.0  0.0 0.00  0.0  0.0  0.0  0.0  0.0   0.0   0.0   0.0  0.00   0.0   0.0  0.50  0.00  0.00  0.00  0.00   0.0  0.00  0.00  0.00  0.00  0.00
 [3,]  0.0  0.0  1.0 0.00  0.0  0.0  0.0  0.0  0.0   0.0   0.0   0.0  0.00   0.0   0.0  0.00  0.00  0.00  0.00  0.00   0.5  0.00  0.00  0.00  0.00  0.00
 [4,]  0.0  0.0  0.0 1.00  0.0  0.0  0.0  0.0  0.0   0.0   0.0   0.0  0.00   0.0   0.5  0.25  0.25  0.00  0.00  0.00   0.0  0.00  0.00  0.00  0.00  0.00
 [5,]  0.0  0.0  0.0 0.00  1.0  0.0  0.0  0.0  0.0   0.0   0.0   0.0  0.00   0.0   0.0  0.00  0.50  0.00  0.00  0.00   0.0  0.00  0.00  0.00  0.00  0.00
 [6,]  0.0  0.0  0.0 0.00  0.0  1.0  0.0  0.0  0.0   0.0   0.0   0.0  0.00   0.0   0.0  0.00  0.00  0.50  0.00  0.00   0.0  0.00  0.00  0.00  0.00  0.00
 [7,]  0.0  0.0  0.0 0.00  0.0  0.0  1.0  0.0  0.0   0.0   0.0   0.0  0.00   0.0   0.0  0.00  0.00  0.00  0.00  0.00   0.0  0.00  0.00  0.00  0.50  0.00
 [8,]  0.0  0.0  0.0 0.00  0.0  0.0  0.0  1.0  0.0   0.0   0.0   0.0  0.00   0.0   0.0  0.00  0.00  0.00  0.00  0.00   0.0  0.50  0.00  0.00  0.00  0.00
 [9,]  0.0  0.0  0.0 0.00  0.0  0.0  0.0  0.0  1.0   0.0   0.0   0.0  0.00   0.0   0.0  0.00  0.00  0.00  0.50  0.50   0.0  0.00  0.00  0.00  0.00  0.00
[10,]  0.0  0.0  0.0 0.00  0.0  0.0  0.0  0.0  0.0   1.0   0.0   0.0  0.00   0.0   0.0  0.00  0.00  0.00  0.00  0.00   0.0  0.00  0.00  0.50  0.00  0.00
[11,]  0.0  0.0  0.0 0.00  0.0  0.0  0.0  0.0  0.0   0.0   1.0   0.0  0.00   0.0   0.0  0.00  0.00  0.00  0.00  0.00   0.0  0.00  0.50  0.00  0.00  0.00
[12,]  0.0  0.0  0.0 0.00  0.0  0.0  0.0  0.0  0.0   0.0   0.0   1.0  0.00   0.0   0.0  0.00  0.00  0.00  0.00  0.00   0.0  0.00  0.00  0.00  0.00  0.50
[13,]  0.0  0.0  0.0 0.00  0.0  0.0  0.0  0.0  0.0   0.0   0.0   0.0  1.00   0.0   0.5  0.25  0.25  0.00  0.00  0.00   0.0  0.00  0.00  0.00  0.00  0.00
[14,]  0.0  0.0  0.0 0.00  0.0  0.0  0.0  0.0  0.0   0.0   0.0   0.0  0.00   1.0   0.0  0.00  0.00  0.50  0.50  0.50   0.0  0.50  0.50  0.50  0.50  0.50
[15,]  0.0  0.0  0.0 0.50  0.0  0.0  0.0  0.0  0.0   0.0   0.0   0.0  0.50   0.0   1.0  0.50  0.50  0.00  0.00  0.00   0.0  0.00  0.00  0.00  0.00  0.00
[16,]  0.0  0.5  0.0 0.25  0.0  0.0  0.0  0.0  0.0   0.0   0.0   0.0  0.25   0.0   0.5  1.00  0.25  0.00  0.00  0.00   0.0  0.00  0.00  0.00  0.00  0.00
[17,]  0.0  0.0  0.0 0.25  0.5  0.0  0.0  0.0  0.0   0.0   0.0   0.0  0.25   0.0   0.5  0.25  1.00  0.00  0.00  0.00   0.0  0.00  0.00  0.00  0.00  0.00
[18,]  0.0  0.0  0.0 0.00  0.0  0.5  0.0  0.0  0.0   0.0   0.0   0.0  0.00   0.5   0.0  0.00  0.00  1.00  0.25  0.25   0.0  0.25  0.25  0.25  0.25  0.25
[19,]  0.0  0.0  0.0 0.00  0.0  0.0  0.0  0.0  0.5   0.0   0.0   0.0  0.00   0.5   0.0  0.00  0.00  0.25  1.00  0.50   0.0  0.25  0.25  0.25  0.25  0.25
[20,]  0.0  0.0  0.0 0.00  0.0  0.0  0.0  0.0  0.5   0.0   0.0   0.0  0.00   0.5   0.0  0.00  0.00  0.25  0.50  1.00   0.0  0.25  0.25  0.25  0.25  0.25
[21,]  0.5  0.0  0.5 0.00  0.0  0.0  0.0  0.0  0.0   0.0   0.0   0.0  0.00   0.0   0.0  0.00  0.00  0.00  0.00  0.00   1.0  0.00  0.00  0.00  0.00  0.00
[22,]  0.0  0.0  0.0 0.00  0.0  0.0  0.0  0.5  0.0   0.0   0.0   0.0  0.00   0.5   0.0  0.00  0.00  0.25  0.25  0.25   0.0  1.00  0.25  0.25  0.25  0.25
[23,]  0.0  0.0  0.0 0.00  0.0  0.0  0.0  0.0  0.0   0.0   0.5   0.0  0.00   0.5   0.0  0.00  0.00  0.25  0.25  0.25   0.0  0.25  1.00  0.25  0.25  0.25
[24,]  0.0  0.0  0.0 0.00  0.0  0.0  0.0  0.0  0.0   0.5   0.0   0.0  0.00   0.5   0.0  0.00  0.00  0.25  0.25  0.25   0.0  0.25  0.25  1.00  0.25  0.25
[25,]  0.0  0.0  0.0 0.00  0.0  0.0  0.5  0.0  0.0   0.0   0.0   0.0  0.00   0.5   0.0  0.00  0.00  0.25  0.25  0.25   0.0  0.25  0.25  0.25  1.00  0.25
[26,]  0.0  0.0  0.0 0.00  0.0  0.0  0.0  0.0  0.0   0.0   0.0   0.5  0.00   0.5   0.0  0.00  0.00  0.25  0.25  0.25   0.0  0.25  0.25  0.25  0.25  1.00
> A22 = A[18:26, 18:26]
> A22
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
 [1,] 1.00 0.25 0.25    0 0.25 0.25 0.25 0.25 0.25
 [2,] 0.25 1.00 0.50    0 0.25 0.25 0.25 0.25 0.25
 [3,] 0.25 0.50 1.00    0 0.25 0.25 0.25 0.25 0.25
 [4,] 0.00 0.00 0.00    1 0.00 0.00 0.00 0.00 0.00
 [5,] 0.25 0.25 0.25    0 1.00 0.25 0.25 0.25 0.25
 [6,] 0.25 0.25 0.25    0 0.25 1.00 0.25 0.25 0.25
 [7,] 0.25 0.25 0.25    0 0.25 0.25 1.00 0.25 0.25
 [8,] 0.25 0.25 0.25    0 0.25 0.25 0.25 1.00 0.25
 [9,] 0.25 0.25 0.25    0 0.25 0.25 0.25 0.25 1.00
> 
> # weighted G
> wt = 0.95
> Gw = wt * G + (1 - wt) * A22
> round(Gw, 3)
        6      7      8      9     10     11     12     13     14
6   0.796  0.130  0.064  0.184 -0.366  0.229 -0.498 -0.101 -0.300
7   0.130  0.730  0.341 -0.742 -0.101  0.196 -0.531  0.163 -0.035
8   0.064  0.341  0.895 -0.511 -0.465  0.725 -0.300 -0.201 -0.399
9   0.184 -0.742 -0.511  2.648 -0.643 -0.048  0.118 -0.676 -0.279
10 -0.366 -0.101 -0.465 -0.643  0.928 -0.598  0.461  0.560  0.362
11  0.229  0.196  0.725 -0.048 -0.598  0.928 -0.432 -0.333 -0.531
12 -0.498 -0.531 -0.300  0.118  0.461 -0.432  0.961  0.130  0.229
13 -0.101  0.163 -0.201 -0.676  0.560 -0.333  0.130  0.565  0.031
14 -0.300 -0.035 -0.399 -0.279  0.362 -0.531  0.229  0.031  1.061
> Gwi = ginv(Gw)
> round(Gwi, 3)
        [,1]   [,2]   [,3]  [,4]   [,5]   [,6]   [,7]   [,8]   [,9]
 [1,]  2.535  0.934  0.907 0.410  1.522  0.065  1.368 -0.812  0.438
 [2,]  0.934  5.634 -1.782 0.567  1.995  1.049  3.148 -3.458 -0.805
 [3,]  0.907 -1.782  6.846 1.825  2.193 -3.183 -1.582  1.542  1.208
 [4,]  0.410  0.567  1.825 1.593  1.227  0.217  0.018  1.324  0.889
 [5,]  1.522  1.995  2.193 1.227  7.496 -0.183 -0.323 -5.482 -0.773
 [6,]  0.065  1.049 -3.183 0.217 -0.183  5.225  1.512  1.685  1.217
 [7,]  1.368  3.148 -1.582 0.018 -0.323  1.512  3.965 -0.903 -0.062
 [8,] -0.812 -3.458  1.542 1.324 -5.482  1.685 -0.903 11.227  3.166
 [9,]  0.438 -0.805  1.208 0.889 -0.773  1.217 -0.062  3.166  2.523
> 
> # inverse matrix of NRM 
> Ai = ainv(pedi) 
> Ai
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26]
 [1,]  1.5  0.0  0.5  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0   0.0     0     0     0     0     0    -1     0     0     0     0     0
 [2,]  0.0  1.5  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0   0.5    -1     0     0     0     0     0     0     0     0     0     0
 [3,]  0.5  0.0  1.5  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0   0.0     0     0     0     0     0    -1     0     0     0     0     0
 [4,]  0.0  0.0  0.0  1.5  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.5   0.0  -1.0     0     0     0     0     0     0     0     0     0     0     0
 [5,]  0.0  0.0  0.0  0.0  1.5  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0   0.5     0    -1     0     0     0     0     0     0     0     0     0
 [6,]  0.0  0.0  0.0  0.0  0.0  1.5  0.0  0.0    0   0.0   0.0   0.0   0.0   0.5   0.0     0     0    -1     0     0     0     0     0     0     0     0
 [7,]  0.0  0.0  0.0  0.0  0.0  0.0  1.5  0.0    0   0.0   0.0   0.0   0.0   0.5   0.0     0     0     0     0     0     0     0     0     0    -1     0
 [8,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.5    0   0.0   0.0   0.0   0.0   0.5   0.0     0     0     0     0     0     0    -1     0     0     0     0
 [9,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    2   0.0   0.0   0.0   0.0   1.0   0.0     0     0     0    -1    -1     0     0     0     0     0     0
[10,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0   1.5   0.0   0.0   0.0   0.5   0.0     0     0     0     0     0     0     0     0    -1     0     0
[11,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0   1.5   0.0   0.0   0.5   0.0     0     0     0     0     0     0     0    -1     0     0     0
[12,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   1.5   0.0   0.5   0.0     0     0     0     0     0     0     0     0     0     0    -1
[13,]  0.0  0.0  0.0  0.5  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   1.5   0.0  -1.0     0     0     0     0     0     0     0     0     0     0     0
[14,]  0.0  0.0  0.0  0.0  0.0  0.5  0.5  0.5    1   0.5   0.5   0.5   0.0   5.0   0.0     0     0    -1    -1    -1     0    -1    -1    -1    -1    -1
[15,]  0.0  0.5  0.0 -1.0  0.5  0.0  0.0  0.0    0   0.0   0.0   0.0  -1.0   0.0   3.0    -1    -1     0     0     0     0     0     0     0     0     0
[16,]  0.0 -1.0  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0  -1.0     2     0     0     0     0     0     0     0     0     0     0
[17,]  0.0  0.0  0.0  0.0 -1.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0  -1.0     0     2     0     0     0     0     0     0     0     0     0
[18,]  0.0  0.0  0.0  0.0  0.0 -1.0  0.0  0.0    0   0.0   0.0   0.0   0.0  -1.0   0.0     0     0     2     0     0     0     0     0     0     0     0
[19,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   -1   0.0   0.0   0.0   0.0  -1.0   0.0     0     0     0     2     0     0     0     0     0     0     0
[20,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   -1   0.0   0.0   0.0   0.0  -1.0   0.0     0     0     0     0     2     0     0     0     0     0     0
[21,] -1.0  0.0 -1.0  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0   0.0     0     0     0     0     0     2     0     0     0     0     0
[22,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0 -1.0    0   0.0   0.0   0.0   0.0  -1.0   0.0     0     0     0     0     0     0     2     0     0     0     0
[23,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0  -1.0   0.0   0.0  -1.0   0.0     0     0     0     0     0     0     0     2     0     0     0
[24,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0  -1.0   0.0   0.0   0.0  -1.0   0.0     0     0     0     0     0     0     0     0     2     0     0
[25,]  0.0  0.0  0.0  0.0  0.0  0.0 -1.0  0.0    0   0.0   0.0   0.0   0.0  -1.0   0.0     0     0     0     0     0     0     0     0     0     2     0
[26,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0  -1.0   0.0  -1.0   0.0     0     0     0     0     0     0     0     0     0     0     2
> 
> # H matrix
> H0 = matrix(c(0), 26, 26)
> H0[18:26, 18:26] = Gwi - ginv(A22)
> H0
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17]      [,18]      [,19]     [,20]      [,21]      [,22]      [,23]       [,24]      [,25]       [,26]
 [1,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  0.0000000  0.0000000  0.000000 0.00000000  0.0000000  0.0000000  0.00000000  0.0000000  0.00000000
 [2,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  0.0000000  0.0000000  0.000000 0.00000000  0.0000000  0.0000000  0.00000000  0.0000000  0.00000000
 [3,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  0.0000000  0.0000000  0.000000 0.00000000  0.0000000  0.0000000  0.00000000  0.0000000  0.00000000
 [4,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  0.0000000  0.0000000  0.000000 0.00000000  0.0000000  0.0000000  0.00000000  0.0000000  0.00000000
 [5,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  0.0000000  0.0000000  0.000000 0.00000000  0.0000000  0.0000000  0.00000000  0.0000000  0.00000000
 [6,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  0.0000000  0.0000000  0.000000 0.00000000  0.0000000  0.0000000  0.00000000  0.0000000  0.00000000
 [7,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  0.0000000  0.0000000  0.000000 0.00000000  0.0000000  0.0000000  0.00000000  0.0000000  0.00000000
 [8,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  0.0000000  0.0000000  0.000000 0.00000000  0.0000000  0.0000000  0.00000000  0.0000000  0.00000000
 [9,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  0.0000000  0.0000000  0.000000 0.00000000  0.0000000  0.0000000  0.00000000  0.0000000  0.00000000
[10,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  0.0000000  0.0000000  0.000000 0.00000000  0.0000000  0.0000000  0.00000000  0.0000000  0.00000000
[11,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  0.0000000  0.0000000  0.000000 0.00000000  0.0000000  0.0000000  0.00000000  0.0000000  0.00000000
[12,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  0.0000000  0.0000000  0.000000 0.00000000  0.0000000  0.0000000  0.00000000  0.0000000  0.00000000
[13,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  0.0000000  0.0000000  0.000000 0.00000000  0.0000000  0.0000000  0.00000000  0.0000000  0.00000000
[14,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  0.0000000  0.0000000  0.000000 0.00000000  0.0000000  0.0000000  0.00000000  0.0000000  0.00000000
[15,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  0.0000000  0.0000000  0.000000 0.00000000  0.0000000  0.0000000  0.00000000  0.0000000  0.00000000
[16,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  0.0000000  0.0000000  0.000000 0.00000000  0.0000000  0.0000000  0.00000000  0.0000000  0.00000000
[17,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  0.0000000  0.0000000  0.000000 0.00000000  0.0000000  0.0000000  0.00000000  0.0000000  0.00000000
[18,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  1.3290432  1.0290474  1.002121 0.40953260  1.6487596  0.1916196  1.49505551 -0.6854842  0.56535657
[19,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  1.0290474  4.2055366 -1.210784 0.56677591  2.0907120  1.1447303  3.24326156 -3.3624887 -0.70927508
[20,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  1.0021206 -1.2107838  5.417277 1.82466493  2.2884902 -3.0879351 -1.48680074  1.6374007  1.30339943
[21,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  0.4095326  0.5667759  1.824665 0.59303100  1.2270063  0.2169856  0.01800482  1.3244139  0.88865550
[22,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  1.6487596  2.0907120  2.288490 1.22700625  6.2894776 -0.0556016 -0.19630790 -5.3550160 -0.64571643
[23,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  0.1916196  1.1447303 -3.087935 0.21698563 -0.0556016  4.0190707  1.63941245  1.8117645  1.34378117
[24,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  1.4950555  3.2432616 -1.486801 0.01800482 -0.1963079  1.6394124  2.75909588 -0.7756486  0.06479731
[25,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0 -0.6854842 -3.3624887  1.637401 1.32441392 -5.3550160  1.8117645 -0.77564863 10.0208796  3.29290470
[26,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0     0     0     0  0.5653566 -0.7092751  1.303399 0.88865550 -0.6457164  1.3437812  0.06479731  3.2929047  1.31668449
> Hi = Ai + H0
> round(Hi, 3)
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17]  [,18]  [,19]  [,20]  [,21]  [,22]  [,23]  [,24]  [,25]  [,26]
 [1,]  1.5  0.0  0.5  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0   0.0     0     0  0.000  0.000  0.000 -1.000  0.000  0.000  0.000  0.000  0.000
 [2,]  0.0  1.5  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0   0.5    -1     0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
 [3,]  0.5  0.0  1.5  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0   0.0     0     0  0.000  0.000  0.000 -1.000  0.000  0.000  0.000  0.000  0.000
 [4,]  0.0  0.0  0.0  1.5  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.5   0.0  -1.0     0     0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
 [5,]  0.0  0.0  0.0  0.0  1.5  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0   0.5     0    -1  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
 [6,]  0.0  0.0  0.0  0.0  0.0  1.5  0.0  0.0    0   0.0   0.0   0.0   0.0   0.5   0.0     0     0 -1.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
 [7,]  0.0  0.0  0.0  0.0  0.0  0.0  1.5  0.0    0   0.0   0.0   0.0   0.0   0.5   0.0     0     0  0.000  0.000  0.000  0.000  0.000  0.000  0.000 -1.000  0.000
 [8,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.5    0   0.0   0.0   0.0   0.0   0.5   0.0     0     0  0.000  0.000  0.000  0.000 -1.000  0.000  0.000  0.000  0.000
 [9,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    2   0.0   0.0   0.0   0.0   1.0   0.0     0     0  0.000 -1.000 -1.000  0.000  0.000  0.000  0.000  0.000  0.000
[10,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0   1.5   0.0   0.0   0.0   0.5   0.0     0     0  0.000  0.000  0.000  0.000  0.000  0.000 -1.000  0.000  0.000
[11,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0   1.5   0.0   0.0   0.5   0.0     0     0  0.000  0.000  0.000  0.000  0.000 -1.000  0.000  0.000  0.000
[12,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   1.5   0.0   0.5   0.0     0     0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 -1.000
[13,]  0.0  0.0  0.0  0.5  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   1.5   0.0  -1.0     0     0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
[14,]  0.0  0.0  0.0  0.0  0.0  0.5  0.5  0.5    1   0.5   0.5   0.5   0.0   5.0   0.0     0     0 -1.000 -1.000 -1.000  0.000 -1.000 -1.000 -1.000 -1.000 -1.000
[15,]  0.0  0.5  0.0 -1.0  0.5  0.0  0.0  0.0    0   0.0   0.0   0.0  -1.0   0.0   3.0    -1    -1  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
[16,]  0.0 -1.0  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0  -1.0     2     0  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
[17,]  0.0  0.0  0.0  0.0 -1.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0  -1.0     0     2  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
[18,]  0.0  0.0  0.0  0.0  0.0 -1.0  0.0  0.0    0   0.0   0.0   0.0   0.0  -1.0   0.0     0     0  3.329  1.029  1.002  0.410  1.649  0.192  1.495 -0.685  0.565
[19,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   -1   0.0   0.0   0.0   0.0  -1.0   0.0     0     0  1.029  6.206 -1.211  0.567  2.091  1.145  3.243 -3.362 -0.709
[20,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0   -1   0.0   0.0   0.0   0.0  -1.0   0.0     0     0  1.002 -1.211  7.417  1.825  2.288 -3.088 -1.487  1.637  1.303
[21,] -1.0  0.0 -1.0  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0   0.0   0.0   0.0   0.0     0     0  0.410  0.567  1.825  2.593  1.227  0.217  0.018  1.324  0.889
[22,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0 -1.0    0   0.0   0.0   0.0   0.0  -1.0   0.0     0     0  1.649  2.091  2.288  1.227  8.289 -0.056 -0.196 -5.355 -0.646
[23,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0  -1.0   0.0   0.0  -1.0   0.0     0     0  0.192  1.145 -3.088  0.217 -0.056  6.019  1.639  1.812  1.344
[24,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0  -1.0   0.0   0.0   0.0  -1.0   0.0     0     0  1.495  3.243 -1.487  0.018 -0.196  1.639  4.759 -0.776  0.065
[25,]  0.0  0.0  0.0  0.0  0.0  0.0 -1.0  0.0    0   0.0   0.0   0.0   0.0  -1.0   0.0     0     0 -0.685 -3.362  1.637  1.324 -5.355  1.812 -0.776 12.021  3.293
[26,]  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0   0.0   0.0  -1.0   0.0  -1.0   0.0     0     0  0.565 -0.709  1.303  0.889 -0.646  1.344  0.065  3.293  3.317
> 
> # weight R
> R = diag(data[1:10, 5])
> R
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]  558    0    0    0    0    0    0    0    0     0
 [2,]    0  722    0    0    0    0    0    0    0     0
 [3,]    0    0  300    0    0    0    0    0    0     0
 [4,]    0    0    0   73    0    0    0    0    0     0
 [5,]    0    0    0    0   52    0    0    0    0     0
 [6,]    0    0    0    0    0   87    0    0    0     0
 [7,]    0    0    0    0    0    0   64    0    0     0
 [8,]    0    0    0    0    0    0    0  103    0     0
 [9,]    0    0    0    0    0    0    0    0   13     0
[10,]    0    0    0    0    0    0    0    0    0   125
> Ri = ginv(R)
> Ri
             [,1]        [,2]        [,3]       [,4]       [,5]       [,6]     [,7]        [,8]       [,9] [,10]
 [1,] 0.001792115 0.000000000 0.000000000 0.00000000 0.00000000 0.00000000 0.000000 0.000000000 0.00000000 0.000
 [2,] 0.000000000 0.001385042 0.000000000 0.00000000 0.00000000 0.00000000 0.000000 0.000000000 0.00000000 0.000
 [3,] 0.000000000 0.000000000 0.003333333 0.00000000 0.00000000 0.00000000 0.000000 0.000000000 0.00000000 0.000
 [4,] 0.000000000 0.000000000 0.000000000 0.01369863 0.00000000 0.00000000 0.000000 0.000000000 0.00000000 0.000
 [5,] 0.000000000 0.000000000 0.000000000 0.00000000 0.01923077 0.00000000 0.000000 0.000000000 0.00000000 0.000
 [6,] 0.000000000 0.000000000 0.000000000 0.00000000 0.00000000 0.01149425 0.000000 0.000000000 0.00000000 0.000
 [7,] 0.000000000 0.000000000 0.000000000 0.00000000 0.00000000 0.00000000 0.015625 0.000000000 0.00000000 0.000
 [8,] 0.000000000 0.000000000 0.000000000 0.00000000 0.00000000 0.00000000 0.000000 0.009708738 0.00000000 0.000
 [9,] 0.000000000 0.000000000 0.000000000 0.00000000 0.00000000 0.00000000 0.000000 0.000000000 0.07692308 0.000
[10,] 0.000000000 0.000000000 0.000000000 0.00000000 0.00000000 0.00000000 0.000000 0.000000000 0.00000000 0.008
> 
> # Left Hand Side
> LHS = rbind(
+   cbind(t(X) %*% Ri %*% X, t(X) %*% Ri %*% W), 
+   cbind(t(W) %*% Ri %*% X, t(W) %*% Ri %*% W + Hi * alpha))
> round(LHS, 3)
       [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]  [,10]  [,11]  [,12]  [,13]  [,14]  [,15]  [,16]  [,17]  [,18]  [,19]   [,20]   [,21]  [,22]   [,23]   [,24]   [,25]   [,26]  [,27]
 [1,] 0.161  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.002  0.001  0.003  0.014  0.019  0.011   0.016   0.010  0.077   0.008   0.000   0.000   0.000  0.000
 [2,] 0.000 10.428  0.000  3.476  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000   0.000   0.000 -6.952   0.000   0.000   0.000   0.000  0.000
 [3,] 0.000  0.000 10.428  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  3.476 -6.952  0.000  0.000   0.000   0.000  0.000   0.000   0.000   0.000   0.000  0.000
 [4,] 0.000  3.476  0.000 10.428  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000   0.000   0.000 -6.952   0.000   0.000   0.000   0.000  0.000
 [5,] 0.000  0.000  0.000  0.000 10.428  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  3.476  0.000 -6.952  0.000  0.000  0.000   0.000   0.000  0.000   0.000   0.000   0.000   0.000  0.000
 [6,] 0.000  0.000  0.000  0.000  0.000 10.428  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  3.476  0.000 -6.952  0.000   0.000   0.000  0.000   0.000   0.000   0.000   0.000  0.000
 [7,] 0.000  0.000  0.000  0.000  0.000  0.000 10.428  0.000  0.000  0.000  0.000  0.000  0.000  0.000  3.476  0.000  0.000  0.000 -6.952   0.000   0.000  0.000   0.000   0.000   0.000   0.000  0.000
 [8,] 0.000  0.000  0.000  0.000  0.000  0.000  0.000 10.428  0.000  0.000  0.000  0.000  0.000  0.000  3.476  0.000  0.000  0.000  0.000   0.000   0.000  0.000   0.000   0.000   0.000  -6.952  0.000
 [9,] 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 10.428  0.000  0.000  0.000  0.000  0.000  3.476  0.000  0.000  0.000  0.000   0.000   0.000  0.000  -6.952   0.000   0.000   0.000  0.000
[10,] 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 13.904  0.000  0.000  0.000  0.000  6.952  0.000  0.000  0.000  0.000  -6.952  -6.952  0.000   0.000   0.000   0.000   0.000  0.000
[11,] 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 10.428  0.000  0.000  0.000  3.476  0.000  0.000  0.000  0.000   0.000   0.000  0.000   0.000   0.000  -6.952   0.000  0.000
[12,] 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 10.428  0.000  0.000  3.476  0.000  0.000  0.000  0.000   0.000   0.000  0.000   0.000  -6.952   0.000   0.000  0.000
[13,] 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 10.428  0.000  3.476  0.000  0.000  0.000  0.000   0.000   0.000  0.000   0.000   0.000   0.000   0.000 -6.952
[14,] 0.002  0.000  0.000  0.000  3.476  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 10.430  0.000 -6.952  0.000  0.000  0.000   0.000   0.000  0.000   0.000   0.000   0.000   0.000  0.000
[15,] 0.001  0.000  0.000  0.000  0.000  0.000  3.476  3.476  3.476  6.952  3.476  3.476  3.476  0.000 34.762  0.000  0.000  0.000 -6.952  -6.952  -6.952  0.000  -6.952  -6.952  -6.952  -6.952 -6.952
[16,] 0.003  0.000  3.476  0.000 -6.952  3.476  0.000  0.000  0.000  0.000  0.000  0.000  0.000 -6.952  0.000 20.860 -6.952 -6.952  0.000   0.000   0.000  0.000   0.000   0.000   0.000   0.000  0.000
[17,] 0.014  0.000 -6.952  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 -6.952 13.918  0.000  0.000   0.000   0.000  0.000   0.000   0.000   0.000   0.000  0.000
[18,] 0.019  0.000  0.000  0.000  0.000 -6.952  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 -6.952  0.000 13.923  0.000   0.000   0.000  0.000   0.000   0.000   0.000   0.000  0.000
[19,] 0.011  0.000  0.000  0.000  0.000  0.000 -6.952  0.000  0.000  0.000  0.000  0.000  0.000  0.000 -6.952  0.000  0.000  0.000 23.155   7.154   6.967  2.847  11.462   1.332  10.394  -4.766  3.930
[20,] 0.016  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 -6.952  0.000  0.000  0.000  0.000 -6.952  0.000  0.000  0.000  7.154  43.157  -8.418  3.940  14.535   7.958  22.548 -23.376 -4.931
[21,] 0.010  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 -6.952  0.000  0.000  0.000  0.000 -6.952  0.000  0.000  0.000  6.967  -8.418  51.576 12.685  15.910 -21.468 -10.336  11.383  9.061
[22,] 0.077 -6.952  0.000 -6.952  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  2.847   3.940  12.685 18.104   8.530   1.509   0.125   9.207  6.178
[23,] 0.008  0.000  0.000  0.000  0.000  0.000  0.000  0.000 -6.952  0.000  0.000  0.000  0.000  0.000 -6.952  0.000  0.000  0.000 11.462  14.535  15.910  8.530  57.638  -0.387  -1.365 -37.229 -4.489
[24,] 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 -6.952  0.000  0.000 -6.952  0.000  0.000  0.000  1.332   7.958 -21.468  1.509  -0.387  41.845  11.397  12.596  9.342
[25,] 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 -6.952  0.000  0.000  0.000 -6.952  0.000  0.000  0.000 10.394  22.548 -10.336  0.125  -1.365  11.397  33.086  -5.392  0.450
[26,] 0.000  0.000  0.000  0.000  0.000  0.000  0.000 -6.952  0.000  0.000  0.000  0.000  0.000  0.000 -6.952  0.000  0.000  0.000 -4.766 -23.376  11.383  9.207 -37.229  12.596  -5.392  83.571 22.893
[27,] 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000 -6.952  0.000 -6.952  0.000  0.000  0.000  3.930  -4.931   9.061  6.178  -4.489   9.342   0.450  22.893 23.058
> 
> # Right Hand Side
> RHS = rbind(t(X) %*% Ri %*% y, 
+             t(W) %*% Ri %*% y)
> round(RHS, 3)
       [,1]
 [1,] 1.351
 [2,] 0.000
 [3,] 0.000
 [4,] 0.000
 [5,] 0.000
 [6,] 0.000
 [7,] 0.000
 [8,] 0.000
 [9,] 0.000
[10,] 0.000
[11,] 0.000
[12,] 0.000
[13,] 0.000
[14,] 0.016
[15,] 0.019
[16,] 0.042
[17,] 0.211
[18,] 0.113
[19,] 0.089
[20,] 0.159
[21,] 0.047
[22,] 0.585
[23,] 0.070
[24,] 0.000
[25,] 0.000
[26,] 0.000
[27,] 0.000
> 
> # Solutions
> 
> # generalized inverse of LHS
> gi_LHS = ginv(LHS)
> round(gi_LHS, 3)
        [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]  [,10]  [,11]  [,12]  [,13]  [,14]  [,15]  [,16]  [,17]  [,18]  [,19]  [,20]  [,21]  [,22]  [,23]  [,24]  [,25]  [,26]  [,27]
 [1,]  6.279 -0.080 -0.006 -0.080 -0.009 -0.009 -0.022  0.019  0.022  0.018 -0.006 -0.011  0.007 -0.011  0.026 -0.019 -0.019 -0.022 -0.020  0.036  0.025 -0.161  0.046 -0.004  0.004  0.042  0.023
 [2,] -0.080  0.201  0.000  0.058  0.000  0.000  0.019 -0.022 -0.020 -0.029  0.016  0.008 -0.003  0.000 -0.031  0.000  0.000  0.000  0.013 -0.052 -0.036  0.187 -0.046 -0.003  0.008 -0.048 -0.020
 [3,] -0.006  0.000  0.144  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.072  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
 [4,] -0.080  0.058  0.000  0.201  0.000  0.000  0.019 -0.022 -0.020 -0.029  0.016  0.008 -0.003  0.000 -0.031  0.000  0.000  0.000  0.013 -0.052 -0.036  0.187 -0.046 -0.003  0.008 -0.048 -0.020
 [5,] -0.009  0.000  0.000  0.000  0.144  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.072  0.036  0.036  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
 [6,] -0.009  0.000  0.000  0.000  0.000  0.144  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.072  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
 [7,] -0.022  0.019  0.000  0.019  0.000  0.000  0.154 -0.004 -0.022  0.015 -0.026  0.022 -0.015  0.000 -0.021  0.000  0.000  0.000  0.077  0.007  0.003  0.038 -0.043  0.022 -0.050 -0.017 -0.033
 [8,]  0.019 -0.022  0.000 -0.022  0.000  0.000 -0.004  0.128  0.031 -0.004  0.008 -0.020  0.000  0.000 -0.004  0.000  0.000  0.000 -0.009  0.010 -0.023 -0.043  0.045 -0.031  0.010  0.046 -0.002
 [9,]  0.022 -0.020  0.000 -0.020  0.000  0.000 -0.022  0.031  0.150 -0.031  0.029 -0.037  0.021  0.000 -0.002  0.000  0.000  0.000 -0.034 -0.016 -0.048 -0.040  0.080 -0.057  0.042  0.046  0.030
[10,]  0.018 -0.029  0.000 -0.029  0.000  0.000  0.015 -0.004 -0.031  0.155 -0.037  0.050 -0.020  0.000 -0.015  0.000  0.000  0.000  0.015  0.068  0.084 -0.057 -0.054  0.067 -0.063 -0.014 -0.038
[11,] -0.006  0.016  0.000  0.016  0.000  0.000 -0.026  0.008  0.029 -0.037  0.161 -0.022  0.016  0.000 -0.015  0.000  0.000  0.000 -0.047 -0.056 -0.032  0.031  0.036 -0.041  0.090  0.005  0.017
[12,] -0.011  0.008  0.000  0.008  0.000  0.000  0.022 -0.020 -0.037  0.050 -0.022  0.162 -0.030  0.000 -0.020  0.000  0.000  0.000  0.023  0.013  0.066  0.016 -0.066  0.089 -0.043 -0.039 -0.056
[13,]  0.007 -0.003  0.000 -0.003  0.000  0.000 -0.015  0.000  0.021 -0.020  0.016 -0.030  0.164  0.000 -0.009  0.000  0.000  0.000 -0.028 -0.009 -0.041 -0.006  0.026 -0.050  0.020 -0.005  0.097
[14,] -0.011  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.144  0.000  0.072  0.036  0.036  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
[15,]  0.026 -0.031  0.000 -0.031  0.000  0.000 -0.021 -0.004 -0.002 -0.015 -0.015 -0.020 -0.009  0.000  0.057  0.000  0.000  0.000 -0.003  0.017  0.009 -0.061  0.025 -0.001  0.006  0.022  0.014
[16,] -0.019  0.000  0.000  0.000  0.072  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.072  0.000  0.144  0.072  0.072  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
[17,] -0.019  0.000  0.072  0.000  0.036  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.036  0.000  0.072  0.144  0.036  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000
[18,] -0.022  0.000  0.000  0.000  0.036  0.072  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.036  0.000  0.072  0.036  0.143  0.000  0.000  0.000  0.001  0.000  0.000  0.000  0.000  0.000
[19,] -0.020  0.013  0.000  0.013  0.000  0.000  0.077 -0.009 -0.034  0.015 -0.047  0.023 -0.028  0.000 -0.003  0.000  0.000  0.000  0.114  0.019  0.009  0.026 -0.052  0.033 -0.072 -0.014 -0.043
[20,]  0.036 -0.052  0.000 -0.052  0.000  0.000  0.007  0.010 -0.016  0.068 -0.056  0.013 -0.009  0.000  0.017  0.000  0.000  0.000  0.019  0.104  0.048 -0.104 -0.015  0.028 -0.076  0.023 -0.005
[21,]  0.025 -0.036  0.000 -0.036  0.000  0.000  0.003 -0.023 -0.048  0.084 -0.032  0.066 -0.041  0.000  0.009  0.000  0.000  0.000  0.009  0.048  0.128 -0.072 -0.067  0.104 -0.043 -0.029 -0.057
[22,] -0.161  0.187  0.000  0.187  0.000  0.000  0.038 -0.043 -0.040 -0.057  0.031  0.016 -0.006  0.000 -0.061  0.000  0.000  0.001  0.026 -0.104 -0.072  0.374 -0.091 -0.007  0.016 -0.095 -0.040
[23,]  0.046 -0.046  0.000 -0.046  0.000  0.000 -0.043  0.045  0.080 -0.054  0.036 -0.066  0.026  0.000  0.025  0.000  0.000  0.000 -0.052 -0.015 -0.067 -0.091  0.133 -0.086  0.066  0.080  0.052
[24,] -0.004 -0.003  0.000 -0.003  0.000  0.000  0.022 -0.031 -0.057  0.067 -0.041  0.089 -0.050  0.000 -0.001  0.000  0.000  0.000  0.033  0.028  0.104 -0.007 -0.086  0.133 -0.062 -0.048 -0.076
[25,]  0.004  0.008  0.000  0.008  0.000  0.000 -0.050  0.010  0.042 -0.063  0.090 -0.043  0.020  0.000  0.006  0.000  0.000  0.000 -0.072 -0.076 -0.043  0.016  0.066 -0.062  0.138  0.019  0.033
[26,]  0.042 -0.048  0.000 -0.048  0.000  0.000 -0.017  0.046  0.046 -0.014  0.005 -0.039 -0.005  0.000  0.022  0.000  0.000  0.000 -0.014  0.023 -0.029 -0.095  0.080 -0.048  0.019  0.081  0.004
[27,]  0.023 -0.020  0.000 -0.020  0.000  0.000 -0.033 -0.002  0.030 -0.038  0.017 -0.056  0.097  0.000  0.014  0.000  0.000  0.000 -0.043 -0.005 -0.057 -0.040  0.052 -0.076  0.033  0.004  0.152
> 
> # solution
> sol = gi_LHS %*% RHS
> 
> # general mean effect
> gme = sol[1]
> round(gme, 3)
[1] 8.39
> 
> # (g)ebv of animals
> gebv = sol[2:27]
> round(gebv, 3)
 [1] -0.012  0.007 -0.012  0.003 -0.003 -0.003  0.004  0.004  0.002 -0.002 -0.003  0.002  0.003  0.004  0.006  0.013 -0.002 -0.002  0.007  0.001 -0.024  0.008 -0.003 -0.001  0.008  0.005
>

결과가 책과 다르다. 책 187쪽에는 5세대의 혈통을 이용하여 구한 A(NRM)가 나오는데 책에서 제시한 혈통을 가지고는 이걸 구할 수가 없다. 본 예제에서 책 187쪽의 혈통을 사용한다고 하였는데 프로그램의 흐름상 사용할 수가 없었다. 이런 이유로 해서 같은 해가 나올 수 없는 것으로 보인다. 본 예제를 위의 참고 포스팅에서는 blupf90으로 푸는데 역시 해가 다르다. 그 이유는 참고 포스팅에서는 EDC를 가중치로 줄 때 그대로 사용하였다. EDC가 아니라 EDC의 역수를 가중치로 주어야 한다. EDC의 역수를 가중치로 주었을 때 동일한 해를 구할 수 있었다.

저작자표시 (새창열림)

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

혈통 그루핑 (0)	2023.05.08
R을 이용한 사회적 관계 모형(Social Interaction Model)의 육종가 구하기 (0)	2021.04.15
GBLUP Model with residual polygenic effects (0)	2020.12.29
SNP-BLUP Model with Polygenic Effects(unweighted analysis) (0)	2020.12.24
Mixed Linear Model(GBLUP model) for Computing DGV (0)	2020.12.21

GBLUP Model with residual polygenic effects

투정이 2020. 12. 29. 21:50

2020. 12. 29. 21:50

# 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를 추가한다. 간단한 설명은 다음 포스팅을 참고한다.

2020/12/24 - [Animal Breeding/R for Genetic Evaluation] - SNP-BLUP Model with Polygenic Effects(unweighted analysis)

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

PREV 이전 1 ···3 4 5 6 7 8 9 ···54 NEXT 다음

Let's Study Animal Breeding

전체 글

R을 이용한 사회적 관계 모형(Social Interaction Model)의 육종가 구하기

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

H 행렬, H 행렬의 역행렬 유도

'Animal Breeding > Genomic Selection' 카테고리의 다른 글

blupf90으로 single-step GBLUP Model 풀기

'Animal Breeding > BLUPF90' 카테고리의 다른 글

Single-step GBLUP(ssGBLUP)

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

GBLUP Model with residual polygenic effects

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

+ Recent posts

티스토리툴바