blupf90으로 fixed regression model의 육종가 구하기

투정이 2014. 3. 14. 01:14

2014. 3. 14. 01:14

blupf90으로 fixed regression model의 육종가 구하기

R. A. Mrode, Linear Models for the prediction of Animal Breeding Values, 2nd Edition. p137 Example 7.1

data

4 1 4 17.0

4 2 38 18.6

4 3 72 24.0

4 4 106 20.0

4 5 140 20.0

4 6 174 15.6

4 7 208 16.0

4 8 242 13.0

4 9 276 8.2

4 10 310 8.0

5 1 4 23.0

5 2 38 21.0

5 3 72 18.0

5 4 106 17.0

5 5 140 16.2

5 6 174 14.0

5 7 208 14.2

5 8 242 13.4

5 9 276 11.8

5 10 310 11.4

6 6 4 10.4

6 7 38 12.3

6 8 72 13.2

6 9 106 11.6

6 10 140 8.4

7 4 4 22.8

7 5 38 22.4

7 6 72 21.4

7 7 106 18.8

7 8 140 18.3

7 9 174 16.2

7 10 208 15.0

8 1 4 22.2

8 2 38 20.0

8 3 72 21.0

8 4 106 23.0

8 5 140 16.8

8 6 174 11.0

8 7 208 13.0

8 8 242 17.0

8 9 276 13.0

8 10 310 12.6

animal, htd(herd-test-day), dim(days-in-milk), test day fat yield

data.txt로 저장

renumf90(blupf90)은 fixed regression model을 위한 다항식(polynomial)을 자동으로 만들어 주지 않는다. 그래서 미리 dim에 대한 다항식을 만들어 주어야 한다.

dim 변수만 뽑는다.

106

140

174

208

242

276

310

106

140

174

208

242

276

310

106

140

106

140

174

208

106

140

174

208

242

276

310

dim.txt로 저장한다.

다음 R 프로그램을 실행시킨다.

# renumf90에선 자동으로 regression model을 위한 polynomial을 생성하지 못함

# renumf90에서 regression model을 적용하기 위한 polynomial 생성 프로그램

# 작업 디렉토리 설정

# 작업 디렉토리를 다음 " " 사이에 복사, 윈도우즈에선 \ 대신에 / 사용

setwd("C:/AB/blupf90/09_fixed_reg")

getwd()

# MASS library 사용

library(MASS)

# read data and check

dim = as.matrix(read.table("dim.txt"))

dim

# maximum

dim_max = max(dim)

dim_max

# minimum

dim_min = min(dim)

dim_min

# make M matrix(polynomials of the standardized DIM values)

dim_1 = c(rep(1,nrow(dim)))

dim_2 = -1 + 2 * (dim - dim_min) / (dim_max - dim_min)

dim_3 = dim_2 ^ 2

dim_4 = dim_3 * dim_2

dim_5 = dim_4 * dim_2

M = cbind(dim_1, dim_2, dim_3, dim_4, dim_5)

# Legendre polynomials 계수를 포함하는 차수 k = 5의 행렬

lambda = matrix(c(

sqrt(.5) , 0.0 , 0.0 , 0.0 , 0.0,

0.0 , sqrt(1.5) , 0.0 , 0.0 , 0.0,

-sqrt(2.5)*.5 , 0.0 , sqrt(2.5)*1.5 , 0.0 , 0.0,

0.0 , -sqrt(3.5)*1.5, 0.0 , sqrt(3.5)*2.5, 0.0,

sqrt(4.5)*3./8., 0.0 , -sqrt(4.5)*30./8., 0.0 , sqrt(4.5)*35./8.),nrow = 5)

lambda

# Legendre Polynomial

dim_poly = M %*% lambda

dim_poly

# export data

write.table(dim_poly, "dim_poly.txt", sep = " ", row.names = FALSE, col.names = FALSE)

위 R 프로그램을 legendre_polynomial.R로 저장

위 프로그램에 주석으로 설명을 추가하였으므로 추가적인 설명은 하지 않음

다음은 R 프로그램 실행화면

“파일 -> 스크립트 열기” 하여 위에서 저장한 legendre_polynomial.R 프로그램을 연다.

dim.txt가 있는 디렉토리를 setwd에 적어준다.

실행할 스크립트 부분을 선택한 후 ctrl + R 을 눌러 실행하면 실행결과가 왼쪽 창에 나타난다.

다음은 R 프로그램 실행결과 생긴 dim_poly.txt의 내용

0.707106781186548 -1.22474487139159 1.58113883008419 -1.87082869338697 2.12132034355964

0.707106781186548 -0.95257934441568 0.644167671515781 -0.0179640615277216 -0.620456293139911

0.707106781186548 -0.680413817439772 -0.0585606974105255 0.757056878668253 -0.775651197104036

0.707106781186548 -0.408248290463863 -0.52704627669473 0.762189467676173 0.0261891400439461

0.707106781186548 -0.136082763487954 -0.761289066336832 0.305389045971261 0.698700390555157

0.707106781186548 0.136082763487954 -0.761289066336832 -0.305389045971261 0.698700390555157

0.707106781186548 0.408248290463863 -0.52704627669473 -0.762189467676173 0.0261891400439465

0.707106781186548 0.680413817439772 -0.0585606974105255 -0.757056878668253 -0.775651197104036

0.707106781186548 0.95257934441568 0.64416767151578 0.0179640615277207 -0.620456293139912

0.707106781186548 1.22474487139159 1.58113883008419 1.87082869338697 2.12132034355964