# Parameter file for program renf90; it is translated to parameter
# file for BLUPF90 family programs.
DATAFILE
fem_snp_data2.txt
TRAITS
6
FIELDS_PASSED TO OUTPUT
WEIGHT(S)
7
RESIDUAL_VARIANCE
245
EFFECT
4 cross alpha
EFFECT
8 cov
EFFECT
9 cov
EFFECT
10 cov
EFFECT
1 cross alpha
RANDOM
animal
FILE
fem_snp_pedi.txt
FILE_POS
1 2 3
PED_DEPTH
0
(CO)VARIANCES
35.241
OPTION solv_method FSPAK
Weight 로 7행을 주었다.
Renumf90 실행화면
Renumf90 실행 로그
RENUMF90 version 1.145
renumf90_fem_snp_w.par
datafile:fem_snp_data2.txt
traits: 6
R
245.0
Processing effect 1 of type cross
item_kind=alpha
Processing effect 2 of type cov
Processing effect 3 of type cov
Processing effect 4 of type cov
Processing effect 5 of type cross
item_kind=alpha
pedigree file name "fem_snp_pedi.txt"
positions of animal, sire, dam, alternate dam, yob, and group 1 2 3 0 0 0 0
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.00179211 1.3571429 -0.3571429 0.2857143
14 0 0 1 722 13.4 0.00138504 0.3571429 -0.3571429 -0.7142857
15 13 4 1 300 12.7 0.00333333 0.3571429 0.6428571 1.2857143
read 8 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 5 of column 1 with 8 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.8875 3.7434 8
8 -0.64286 1.3571 -0.17857E-01 0.74402 8
9 -0.35714 0.64286 0.14286 0.53452 8
10 -0.71429 1.2857 0.28571 0.75593 8
Correlation matrix
6 8 9 10
6 1.00 0.12 -0.60 0.35
8 0.12 1.00 -0.18 -0.25
9 -0.60 -0.18 1.00 0.00
10 0.35 -0.25 0.00 1.00
Counts of nonzero values (order as above)
8 8 8 8
8 8 8 8
8 8 8 8
8 8 8 8
random effect 5
type:animal
opened output pedigree file "renadd05.ped"
read 26 pedigree records
loaded 18 parent(s) in round 0
Pedigree checks
Number of animals with records = 8
Number of parents without records = 18
Total number of animals = 26
Wrote parameter file "renf90.par"
Wrote renumbered data "renf90.dat" 8 records
Renumf90 실행 결과 생성 파일
renf90.tables
Effect group 1 of column 1 with 1 levels, effect # 1
Value # consecutive number
1 8 1
renf90.par
BLUPF90 ver. 1.68
Parameter file: renf90.par
Data file: renf90.dat
Number of Traits 1
Number of Effects 5
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 covariable 4 1
3 covariable 5 1
4 covariable 6 1
5 cross-classified 7 26
Residual (co)variance Matrix
245.00
Random Effect(s) 5
Type of Random Effect: additive animal
Pedigree File: renadd05.ped
trait effect (CO)VARIANCES
1 5 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
Data record length = 7
# equations = 30
G
35.241
read 8 records in 4.6875000E-02 s, 50
nonzeroes
read 26 additive pedigrees
finished peds in 4.6875000E-02 s, 103 nonzeroes
solutions stored in file: "solutions"
# Parameter file for program renf90; it is translated to parameter
# file for BLUPF90 family programs.
DATAFILE
fem_snp_data2.txt
TRAITS
6
FIELDS_PASSED TO OUTPUT
WEIGHT(S)
RESIDUAL_VARIANCE
245
EFFECT
4 cross alpha
EFFECT
8 cov
EFFECT
9 cov
EFFECT
10 cov
EFFECT
1 cross alpha
RANDOM
animal
FILE
fem_snp_pedi.txt
FILE_POS
1 2 3
PED_DEPTH
0
(CO)VARIANCES
35.241
OPTION solv_method FSPAK
실행
명령창에서 다음과 같은 명령어로 renumf90을 실행한다.
renumf90 renumf90_fem_snp_uw.par | tee renumf90_fem_snp_uw_01.log
Renumf90 실행 로그
RENUMF90 version 1.145
renumf90_fem_snp_uw.par
datafile:fem_snp_data2.txt
traits: 6
R
245.0
Processing effect 1 of type cross
item_kind=alpha
Processing effect 2 of type cov
Processing effect 3 of type cov
Processing effect 4 of type cov
Processing effect 5 of type cross
item_kind=alpha
pedigree file name "fem_snp_pedi.txt"
positions of animal, sire, dam, alternate dam, yob, and group 1 2 3 0 0 0 0
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.00179211 1.3571429 -0.3571429 0.2857143
14 0 0 1 722 13.4 0.00138504 0.3571429 -0.3571429 -0.7142857
15 13 4 1 300 12.7 0.00333333 0.3571429 0.6428571 1.2857143
read 8 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 5 of column 1 with 8 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.8875 3.7434 8
8 -0.64286 1.3571 -0.17857E-01 0.74402 8
9 -0.35714 0.64286 0.14286 0.53452 8
10 -0.71429 1.2857 0.28571 0.75593 8
Correlation matrix
6 8 9 10
6 1.00 0.12 -0.60 0.35
8 0.12 1.00 -0.18 -0.25
9 -0.60 -0.18 1.00 0.00
10 0.35 -0.25 0.00 1.00
Counts of nonzero values (order as above)
8 8 8 8
8 8 8 8
8 8 8 8
8 8 8 8
random effect 5
type:animal
opened output pedigree file "renadd05.ped"
read 26 pedigree records
loaded 18 parent(s) in round 0
Pedigree checks
Number of animals with records = 8
Number of parents without records = 18
Total number of animals = 26
Wrote parameter file "renf90.par"
Wrote renumbered data "renf90.dat" 8 records
Renumf90 실형 결과 파일
renf90.tables
Effect group 1 of column 1 with 1 levels, effect # 1
Value # consecutive number
1 8 1
blupf90 renf90.par | tee blupf90_fem_snp_uw_01.log
실행 화면
blupf90 실행 로그
renf90.par
BLUPF90 ver. 1.68
Parameter file: renf90.par
Data file: renf90.dat
Number of Traits 1
Number of Effects 5
Position of Observations 1
Position of Weight (1) 0
Value of Missing Trait/Observation 0
EFFECTS
# type position (2) levels [positions for nested]
1 cross-classified 2 1
2 covariable 3 1
3 covariable 4 1
4 covariable 5 1
5 cross-classified 6 26
Residual (co)variance Matrix
245.00
Random Effect(s) 5
Type of Random Effect: additive animal
Pedigree File: renadd05.ped
trait effect (CO)VARIANCES
1 5 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
Data record length = 6
# equations = 30
G
35.241
read 8 records in 6.2500000E-02 s, 50
nonzeroes
read 26 additive pedigrees
finished peds in 6.2500000E-02 s, 103 nonzeroes
solutions stored in file: "solutions"
# Linear Models for the Prediction of Animal Breeding Values, 3rd Edition.
# Raphael Mrode
# Example 11.1 p180
SNP effect를 고정효과로 모형에 넣어 SNP 추정하고 개체의 genotype을 이용하여 개체의 DGV(direct genomic breeding value)를 추정할 수 있다. SNP effect로 개체의 모든 유전적 효과를 설명할 수 없으므로 polygenic effect를 추가하여 분석한다. 그래서 분석모형은 y = mean + snp effect + polygenic effect + e 가 된다.
SNP effect의 design matrix를 만드는 방법은 다음과 같다. SNP gentype을 먼저 0, 1, 2로 코드화 한다. SNP 별로 평균을 구하는데 각 개체별로 두 개의 SNP가 동원된 것으므로 SNP 별 합을 개체의 2배수로 나누어 평균을 구한다. 코드화한 genotype에서 각 SNP의 평균의 2배를 빼 SNP effect의 design matrix를 구한다. 이 개체별 design matrix를 이용하여 나중에 SNP effect 추정 후 DGV를 추정한다.
참고 : Linear Models for the Prediction of Animal Breeding Values, 3rd Edtion. 9.3 Random Regression Model
임의 회귀 모형은 해를 구하고 나서도 해야할 일이 많다. 개체의 육종가 그래프를 그려 개체별로 비유기 동안 육종가가 어떻게 변하는지 살펴볼 수도 있다. 사실 모든 개체에 대해서 이렇게 그래프를 그려도 비교하기가 어렵기 때문에 사실상 할 일은 아니지만 그래도 개체의 일일 육종가를 구하여 그래프를 그리는 연습을 해 보자. 왜? 책에 그래프가 있으므로. 그리고 6일부터 310일까지의 일일 육종가를 구하면 305일 육종가가 된다. 아마도 이걸 비교하여 개체를 선발할 것이다. 어떤 프로그램을 사용하든 임의 회귀 모형의 해를 구했다고 가정하고 시작하자. 먼저 비유기 동안 분산의 변화부터 살펴보고나서 개체 육종가 그래프와 305일 육종가를 계산하자.
주의해야 할 것이 하나 있다. 예를 들어 유전평가를 할 때 DIM(days in milk, 분만후 비유일수)을 4에서 310을 했을 경우 즉 DIM 4에서 310까지의 르장드르 다항식을 사용하였을 경우 뒤의 모든 분석에서 그 다항식을 사용해야 한다. 예를 들어 305일 유지량을 계산하려면 DIM 6에서 DIM 310까지의 매일 매일 육종가를 더하는데 이때 DIM 6에서 310까지의 르장드르 다항식을 새로 만들어 계산하면 안된다. 유전 평가에서 사용했던 DIM 4에서 310의 르장드르 다항식 중 6에서 310까지만 추출해서 써야 한다. 책에는 설명이 안 되어 있고 해 보니 깨닫게 되는 중요한 사항이다.
여기서 르장드르 다항식(Legendre Polynomials)을 구하는 것이 중요한데 그것은 인터넷에서 복사한 것임을 밝힌다. 다음 사이트를 참조한다.
# plotting the genetic and permanent environmental variance for DIM
# plotting the animals breeding values for DIM
## Given time points covariate and order of fit for Legendre polynomials, return matrix 'M' containing the polynomials of standardized time. 'M' is order t (number of time points) by k (order of Legendre polynomials).
## Arguments
## t: a vector of time, age or days
## n: order of polynomials
## tmax: max time (optional)
## tmin: min time (optional)
## Literature: Mrode, R.A. 2005. Linear Models for the Prediction of Animal Breeding Values. CAB International, Oxon, UK.
## Author: Gota Morota <morota at wisc dot edu>
## Created: 31-Mar-2010
## Last-Modified: 2-Apr-2010
## License: GPLv3 or later
`stdtime` <-
function(t, n, tmax, tmin){
if(missing(tmax)) {
tmax <- t[which.max(t)]
}
if(missing(tmin)) {
tmin <- t[which.min(t)]
}
N <- n+1
M <- matrix(0, nrow=length(t), ncol=N)
a <- -1 + 2*(t-tmin)/(tmax - tmin)
M[,1] <- 1
for (i in 2:N){
M[,i] <- a^(i-1)
}
return(M)
}
## Return coefficient matrix (lambda) of n-th order Legendre polynomials
## Arguments
## n: order of polynomials
## gengler: logical value. If TRUE, Genlger's scaling (1999) will be applied. If not specified, TRUE is assumed.
## Literatures
## Mrode, R.A. 2005. Linear Models for the Prediction of Animal Breeding Values. CAB International, Oxon, UK.
## Gengler, N. et. al. 1999. Estimation of (Co)variance Function Coefficients for Test Day Yield with a Expectation-Maximization Restricted Maximum Likelihood Algorithm. Journal of Dairy Science. 82
## Author: Gota Morota <morota at wisc dot edu>
## Create: 31-Mar-2010
## Last-Modified: 2-Apr-2010
## License: GPLv3 or later
`legendre` <-
function(n, gengler){
if (nargs()==1){
gengler <- TRUE
}
if (gengler != TRUE & gengler != FALSE){
gengler=TRUE
}
N <- n+1
L <- matrix(0,nrow=N, ncol=N)
for(i in (1:N)){
if(i==1){
L[i,i] <- 1
}
else if(i==2){
L[i,i] <- 1
}
else {
tmp <- L[i-1,]
tmp2 <- as.numeric()
tmp2 <- c(0,tmp[1:(N-1)])
L[i,] <- (1/(i-2+1))*( (2*(i-2) + 1)*tmp2 -(i-2)*L[i-2,] )
}
}
# Normalize
for (j in (1:N)){
L[j,] <- (sqrt( (2*(j-1)+1)/2) )*L[j,]
}
# Gengler (1999)
if (gengler==TRUE){
L <- sqrt(2)*L
}
return(L)
}
## example
## DIM <- c(4,38,72,106,140,174,208,242,276,310)
## order <- 4
## M <- stdtime(DIM, order)
## Lambda <- legendre(order, gengler=FALSE)
## Phi <- M%*%t(Lambda)
## Phi
# days in milk
dim = c(4:310)
head(dim)
tail(dim)
# order of animal regression coefficients
order = 2
# Mu : the matrix containing the polynomials of the standardized DIM values
M <- stdtime(dim, order)
head(M)
tail(M)
## Lambda : a matrix of order k containing the coefficients of Legendre polynomials.
Lambda <- legendre(order, gengler=FALSE)
head(Lambda)
tail(Lambda)
## Phi : Legendre polynomials of days in milk
Phi <- M %*% t(Lambda)
head(Phi)
tail(Phi)
# genetic variance and covariance
G = matrix(c( 3.297, 0.594, -1.381,
0.594, 0.921, -0.289,
-1.321, -0.289, 1.005 ), nrow = 3, byrow = TRUE)
G
P = matrix(c( 6.872, -0.254, -1.101,
-0.254, 3.171, 0.167,
-1.101, 0.167, 2.457 ), nrow = 3, byrow = TRUE)
P
# genetic variance for DIM i (vii)
fdim = 3
vii = matrix(rep(0, 3 * 307) , ncol = 3)
head(vii)
for (i in 1 : 307){
vii[i, 1] = fdim + i
vii[i, 2] = Phi[i,] %*% G %*% matrix(Phi[i,])
vii[i, 3] = Phi[i,] %*% P %*% matrix(Phi[i,])
}
head(vii)
tail(vii)
plot(vii[,1], vii[,2], ylim = c(0,12))
par(new = TRUE)
plot(vii[,1], vii[,3], ylim = c(0,12))
# breeding value of each animal for DIM i
# regression coefficients for animal effects
bvreg = matrix(c(
1, -0.0583, 0.0552, -0.0442,
2, -0.0728, -0.0305, -0.0244,
3, 0.1312, -0.0247, 0.0685,
4, 0.3445, 0.0063, -0.3164,
5, -0.4538, -0.052 , 0.2798,
6, -0.5486, 0.073 , 0.1946,
7, 0.8518, -0.0095, -0.3131,
8, 0.2208, 0.0127, -0.0174
), nrow = 8, byrow = TRUE)
bvreg
bvdim = matrix(rep(0, 9 * 307), ncol = 9)
head(bvdim)
for(i in 1 : 307){
bvdim[i, 1] = fdim + i
for(j in 1 : 8){
bvdim[i, j + 1] = Phi[i,] %*% matrix(bvreg[j, 2:4])
}
}
head(bvdim)
plot(bvdim[,1], bvdim[,3], ylim = c(-0.15,0.25))
par(new = TRUE)
plot(bvdim[,1], bvdim[,4], ylim = c(-0.15,0.25))
par(new = TRUE)
plot(bvdim[,1], bvdim[,9], ylim = c(-0.15,0.25))
# 305d breeding value
Phi2 = Phi[3 : 307,]
dim(Phi2)
t = colSums(Phi2)
t
bv_animal = rep(0, 8)
bv_animal
for(i in 1:8){
bv_animal[i] = t %*% matrix(bvreg[i, 2:4])
}
bv_animal
실행 결과
> # plotting the genetic and permanent environmental variance for DIM
> # plotting the animals breeding values for DIM
>
>
> ## Given time points covariate and order of fit for Legendre polynomials, return matrix 'M' containing the polynomials of standardized time. 'M' is order t (number of time points) by k (order of Legendre polynomials).
>
> ## Arguments
> ## t: a vector of time, age or days
> ## n: order of polynomials
> ## tmax: max time (optional)
> ## tmin: min time (optional)
>
> ## Literature: Mrode, R.A. 2005. Linear Models for the Prediction of Animal Breeding Values. CAB International, Oxon, UK.
>
> ## Author: Gota Morota <morota at wisc dot edu>
> ## Created: 31-Mar-2010
> ## Last-Modified: 2-Apr-2010
> ## License: GPLv3 or later
>
> `stdtime` <-
+ function(t, n, tmax, tmin){
+ if(missing(tmax)) {
+ tmax <- t[which.max(t)]
+ }
+ if(missing(tmin)) {
+ tmin <- t[which.min(t)]
+ }
+
+ N <- n+1
+ M <- matrix(0, nrow=length(t), ncol=N)
+ a <- -1 + 2*(t-tmin)/(tmax - tmin)
+ M[,1] <- 1
+
+ for (i in 2:N){
+ M[,i] <- a^(i-1)
+ }
+
+ return(M)
+ }
>
> ## Return coefficient matrix (lambda) of n-th order Legendre polynomials
>
> ## Arguments
> ## n: order of polynomials
> ## gengler: logical value. If TRUE, Genlger's scaling (1999) will be applied. If not specified, TRUE is assumed.
>
> ## Literatures
> ## Mrode, R.A. 2005. Linear Models for the Prediction of Animal Breeding Values. CAB International, Oxon, UK.
> ## Gengler, N. et. al. 1999. Estimation of (Co)variance Function Coefficients for Test Day Yield with a Expectation-Maximization Restricted Maximum Likelihood Algorithm. Journal of Dairy Science. 82
>
> ## Author: Gota Morota <morota at wisc dot edu>
> ## Create: 31-Mar-2010
> ## Last-Modified: 2-Apr-2010
> ## License: GPLv3 or later
>
> `legendre` <-
+ function(n, gengler){
+
+ if (nargs()==1){
+ gengler <- TRUE
+ }
+
+ if (gengler != TRUE & gengler != FALSE){
+ gengler=TRUE
+ }
+
+ N <- n+1
+ L <- matrix(0,nrow=N, ncol=N)
+
+ for(i in (1:N)){
+ if(i==1){
+ L[i,i] <- 1
+ }
+ else if(i==2){
+ L[i,i] <- 1
+ }
+ else {
+ tmp <- L[i-1,]
+ tmp2 <- as.numeric()
+ tmp2 <- c(0,tmp[1:(N-1)])
+ L[i,] <- (1/(i-2+1))*( (2*(i-2) + 1)*tmp2 -(i-2)*L[i-2,] )
+ }
+ }
+
+ # Normalize
+ for (j in (1:N)){
+ L[j,] <- (sqrt( (2*(j-1)+1)/2) )*L[j,]
+ }
+
+
+ # Gengler (1999)
+ if (gengler==TRUE){
+ L <- sqrt(2)*L
+ }
+
+ return(L)
+
+ }
>
> ## example
> ## DIM <- c(4,38,72,106,140,174,208,242,276,310)
> ## order <- 4
> ## M <- stdtime(DIM, order)
> ## Lambda <- legendre(order, gengler=FALSE)
> ## Phi <- M%*%t(Lambda)
> ## Phi
>
>
> # days in milk
> dim = c(4:310)
> head(dim)
[1] 4 5 6 7 8 9
> tail(dim)
[1] 305 306 307 308 309 310
>
> # order of animal regression coefficients
> order = 2
>
> # Mu : the matrix containing the polynomials of the standardized DIM values
> M <- stdtime(dim, order)
> head(M)
[,1] [,2] [,3]
[1,] 1 -1.0000000 1.0000000
[2,] 1 -0.9934641 0.9869708
[3,] 1 -0.9869281 0.9740271
[4,] 1 -0.9803922 0.9611688
[5,] 1 -0.9738562 0.9483959
[6,] 1 -0.9673203 0.9357085
> tail(M)
[,1] [,2] [,3]
[302,] 1 0.9673203 0.9357085
[303,] 1 0.9738562 0.9483959
[304,] 1 0.9803922 0.9611688
[305,] 1 0.9869281 0.9740271
[306,] 1 0.9934641 0.9869708
[307,] 1 1.0000000 1.0000000
>
> ## Lambda : a matrix of order k containing the coefficients of Legendre polynomials.
> Lambda <- legendre(order, gengler=FALSE)
> head(Lambda)
[,1] [,2] [,3]
[1,] 0.7071068 0.000000 0.000000
[2,] 0.0000000 1.224745 0.000000
[3,] -0.7905694 0.000000 2.371708
> tail(Lambda)
[,1] [,2] [,3]
[1,] 0.7071068 0.000000 0.000000
[2,] 0.0000000 1.224745 0.000000
[3,] -0.7905694 0.000000 2.371708
>
> ## Phi : Legendre polynomials of days in milk
> Phi <- M %*% t(Lambda)
> head(Phi)
[,1] [,2] [,3]
[1,] 0.7071068 -1.224745 1.581139
[2,] 0.7071068 -1.216740 1.550237
[3,] 0.7071068 -1.208735 1.519539
[4,] 0.7071068 -1.200730 1.489043
[5,] 0.7071068 -1.192725 1.458749
[6,] 0.7071068 -1.184721 1.428658
> tail(Phi)
[,1] [,2] [,3]
[302,] 0.7071068 1.184721 1.428658
[303,] 0.7071068 1.192725 1.458749
[304,] 0.7071068 1.200730 1.489043
[305,] 0.7071068 1.208735 1.519539
[306,] 0.7071068 1.216740 1.550237
[307,] 0.7071068 1.224745 1.581139
>
> # genetic variance and covariance
> G = matrix(c( 3.297, 0.594, -1.381,
+ 0.594, 0.921, -0.289,
+ -1.321, -0.289, 1.005 ), nrow = 3, byrow = TRUE)
> G
[,1] [,2] [,3]
[1,] 3.297 0.594 -1.381
[2,] 0.594 0.921 -0.289
[3,] -1.321 -0.289 1.005
>
> P = matrix(c( 6.872, -0.254, -1.101,
+ -0.254, 3.171, 0.167,
+ -1.101, 0.167, 2.457 ), nrow = 3, byrow = TRUE)
> P
[,1] [,2] [,3]
[1,] 6.872 -0.254 -1.101
[2,] -0.254 3.171 0.167
[3,] -1.101 0.167 2.457
>
> # genetic variance for DIM i (vii)
> fdim = 3
> vii = matrix(rep(0, 3 * 307) , ncol = 3)
> head(vii)
[,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
>
> for (i in 1 : 307){
+ vii[i, 1] = fdim + i
+ vii[i, 2] = Phi[i,] %*% G %*% matrix(Phi[i,])
+ vii[i, 3] = Phi[i,] %*% P %*% matrix(Phi[i,])
+ }
> head(vii)
[,1] [,2] [,3]
[1,] 4 2.612026 11.66624
[2,] 5 2.533496 11.42854
[3,] 6 2.457661 11.19690
[4,] 7 2.384484 10.97121
[5,] 8 2.313922 10.75139
[6,] 9 2.245937 10.53735
> tail(vii)
[,1] [,2] [,3]
[302,] 305 2.279769 10.81685
[303,] 306 2.306494 11.05675
[304,] 307 2.334957 11.30292
[305,] 308 2.365192 11.55545
[306,] 309 2.397234 11.81442
[307,] 310 2.431118 12.07994
>
> plot(vii[,1], vii[,2], ylim = c(0,12))
> par(new = TRUE)
> plot(vii[,1], vii[,3], ylim = c(0,12))
>
>
> # breeding value of each animal for DIM i
> # regression coefficients for animal effects
> bvreg = matrix(c(
+ 1, -0.0583, 0.0552, -0.0442,
+ 2, -0.0728, -0.0305, -0.0244,
+ 3, 0.1312, -0.0247, 0.0685,
+ 4, 0.3445, 0.0063, -0.3164,
+ 5, -0.4538, -0.052 , 0.2798,
+ 6, -0.5486, 0.073 , 0.1946,
+ 7, 0.8518, -0.0095, -0.3131,
+ 8, 0.2208, 0.0127, -0.0174
+ ), nrow = 8, byrow = TRUE)
> bvreg
[,1] [,2] [,3] [,4]
[1,] 1 -0.0583 0.0552 -0.0442
[2,] 2 -0.0728 -0.0305 -0.0244
[3,] 3 0.1312 -0.0247 0.0685
[4,] 4 0.3445 0.0063 -0.3164
[5,] 5 -0.4538 -0.0520 0.2798
[6,] 6 -0.5486 0.0730 0.1946
[7,] 7 0.8518 -0.0095 -0.3131
[8,] 8 0.2208 0.0127 -0.0174
>
> bvdim = matrix(rep(0, 9 * 307), ncol = 9)
> head(bvdim)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 0 0 0 0 0 0 0 0 0
[2,] 0 0 0 0 0 0 0 0 0
[3,] 0 0 0 0 0 0 0 0 0
[4,] 0 0 0 0 0 0 0 0 0
[5,] 0 0 0 0 0 0 0 0 0
[6,] 0 0 0 0 0 0 0 0 0
>
> for(i in 1 : 307){
+ bvdim[i, 1] = fdim + i
+ for(j in 1 : 8){
+ bvdim[i, j + 1] = Phi[i,] %*% matrix(bvreg[j, 2:4])
+ }
+ }
> head(bvdim)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 4 -0.1787166 -0.05270244 0.2313316 -0.2643899 0.1852043 -0.1696355 0.1188941 0.1130631
[2,] 5 -0.1769089 -0.05219260 0.2290172 -0.2545623 0.1761419 -0.1750646 0.1284932 0.1137024
[3,] 6 -0.1751101 -0.05168770 0.2267166 -0.2447988 0.1671361 -0.1804542 0.1380290 0.1143383
[4,] 7 -0.1733203 -0.05118774 0.2244299 -0.2350994 0.1581870 -0.1858044 0.1475013 0.1149706
[5,] 8 -0.1715395 -0.05069272 0.2221570 -0.2254641 0.1492946 -0.1911152 0.1569101 0.1155993
[6,] 9 -0.1697676 -0.05020266 0.2198981 -0.2158929 0.1404590 -0.1963865 0.1662555 0.1162246
>
> plot(bvdim[,1], bvdim[,3], ylim = c(-0.15,0.25))
> par(new = TRUE)
> plot(bvdim[,1], bvdim[,4], ylim = c(-0.15,0.25))
> par(new = TRUE)
> plot(bvdim[,1], bvdim[,9], ylim = c(-0.15,0.25))
>
> # 305d breeding value
> Phi2 = Phi[3 : 307,]
> dim(Phi2)
[1] 305 3
> t = colSums(Phi2)
> t
[1] 215.667568 2.441485 -1.545070
>
> bv_animal = rep(0, 8)
> bv_animal
[1] 0 0 0 0 0 0 0 0
>
> for(i in 1:8){
+ bv_animal[i] = t %*% matrix(bvreg[i, 2:4])
+ }
>
> bv_animal
[1] -12.37036 -15.73736 28.12944 74.80172 -98.42921 -118.43767 184.16620 47.67729
>
