# 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
혈통
1 0 0
2 0 0
3 0 0
4 0 0
5 0 0
6 0 0
7 1 4
8 1 4
9 2 5
10 1 4
11 2 5
12 3 6
13 2 5
14 3 6
15 3 6
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 |