Let's Study Animal Breeding

log likelihood(이것에 대한 부연 설명 생략)를 이용하여 분산 성분을 추정하는 세 가지 방법 1) 미분 없이 derivative free 2) 1차 미분 EM-REMl 3) 2차 미분(AI-REML)

AI-REML 설명을 위하여 observed information matrix, expected information matrix 및 average information matrix를 이전 포스트에서 설명하였다. 여기서는 average information matrix(2차 미분)를 이용하여 분산 성분을 추정하는 과정과 그 과정에서 Newton-Raphson 이라는 방법이 동원되는데 그것에 대해서 간략히 설명한다.

먼저 Newton-Raphson 방법은 해를 풀어서 구하기 힘들 때 점근적으로 (여러 라운드에 걸쳐서) 해를 구하는 방법이다. 현재의 해(또는 초기값)를 Xn이라고 하고 다음 번 해를 Xn+1이라고 할 때 새로운 해를 구할 때(또는 갱신할 때) Newton-Raphson 방법을 이용한다는 뜻이다. 해가 여러 개인 다차원 문제는 그림으로 이해하기 어려우므로 해가 1개인 방정식(f(x))을 2차원적으로 이해해 보자. 초기 해(Xn 또는 X0)가 있을 때, 방정식의 y 즉 f(Xn)을 구한다. 점(Xn, f(Xn))은 방정식의 위 점이 된다. 이 점에서 이 방정식의 접선의 식을 구해보자. 그러러면 기울기가 필요한데 미분을 한 식에 Xn을 대입하면 된다. 미분한 식을 f'(X)라 하면 접선의 기울기는 f'(Xn)이 된다. 그러면 접선의 식은 y - f(Xn) = f'(Xn)(X - Xn)이 된다. 이 접선이 y축과 만나는 점을 (Xn+1, 0)이라 하자. 이 점을 접선에 대입하면 0 - f(Xn) = f'(Xn)(Xn+1 - Xn)이 된다. 이것을 정리하면 Xn+1 = Xn - f(Xn)/f'(Xn) 이 된다. 이게 Newton-Raphson 식인데 설명을 하자면 현재값 Xn을 갱신하는 방법이 (현재 방정식 값)/(1차 미분 방정식 값)을 빼주는 방식이다.

위의 그림을 보면 만일 그래프가 뾰족하면 즉 접선의 기울기 절대값이 크면 갱신한 Xn+1이 그리 크게 변하지 않을 것이고, 그래프가 좀 완만하다면 즉 접선의 기울기 절대값이 작다면 갱신한 Xn+1이 크게 변한다는 것을 쉽게 알 수 있다. 그래서 AI-REML이 수렴값 근처에서 즉 점점 평평해 질 때 빠른 속으로 해에 더 가까워 진다는 것이고, EM-REML이 수렴의 방향만을 고려한다면 AI-REML은 방향뿐만 아니라 일종의 폭도 고려한다는 얘기다.

이제 다시 분산 성분 추정하는 문제로 돌아봐 보자. 우리가 1차 미분한 식을 0으로 놓고 해를 구한다고 했는데, 이렇게 된 해에 약간의 보정을 가하는 것이다. 즉 Newton-Raphson 방법을 동원하여 약간의 보정을 가하는 것이다. 그런데 위에서 Netwon-Raphson 방법을 이용하려면 미분을 해야하는데, 우리는 이미 1차 미분을 한 식을 이용하니, 여기에 한 번 더 미분 즉 2차 미분을 한 것이다. 위의 Newton-Raphson 식에서 f(Xn)이 1차 미분 log likelihood이고, f'(Xn)이 2차 미분 log likilihood, 즉 average information matrix이다. 그런데 위에서는 f(Xn)이나 f'(Xn+1)이 스칼라이므로 나누기를 했지만 우리가 다루는 분산 성분 추정에서는 스칼라가 아니라 vector이고 matrix가 된다. 다음과 같이 나누기가 아니라 벡터와 역행렬을 이용한 표기로 바뀐다.

분산 성분 갱신 공식은 위와 같은데 세타는 t, t + 1번째 분산 성분 벡터, M은 average information matrix의 역행렬, d는 경사 벡터(gradient vector) 또는 스코어 함수(score function)으로 1차 미분한 값이다. d가 방향이라면 M은 방향의 크기라고 할 수 있다.

마땅한 예제가 없어 예를 들어 보일 수는 없으나 AI-REML  방법은 다음과 같이 진행된다.

1) 초기값 설정

2) 육종가 추정(일종의 expectation)

3) 분산 성분 추정(일종의 maximization, 즉 여기까지는 EM-REML과 동일)

4) average information matrix 계산

5) Newton-Raphson 방법에 따라 분산 성분 업데이트

6) 수렴 여부 판단

이게 각 round 마다 벌어지는 일이다. 수렴이 되었으면 끝이지만 수렴이 안 되었으면 다시 육종가 추정스텝으로 간다. 수렴이 될 때까지 반복한다. 4번과 5번을 계산하느라 한 round의 시간은 길어지지만 전체적인 round가 줄어서 추정 시간이 준다는게 Gilmour 논문의 주장이다. 하나마나한 이야기지만 분석하는 형질의 수나, 모형의 복잡성에 달라질 것이다. 

observed information matrix(OI)와 expected information matrix(EI)를 이용하여 average information matrix를 만들 차례이다. 이름에 average가 들어가 있으므로 OI와 EI의 평균일 것이라 생각할 수 있다. 그러나 그렇지 않다. Gilmour의 Average Information REML An Efficient Algorithm for Variance Parameter Estimation in Linear Mixed Models 논문에는 simplified average 라고 나와 있다. 이게 무슨 뜻이냐. trace는 계산하기 힘들므로 trace를 quadratic form(y'Qy)으로 대체(또는 근사?)한 후 평균을 구했다는 뜻이다. 그렇게 해야 계산도 빨라지고, positive semidefinite이 보장되어 계산이 된다는 뜻이다. 즉 빠르고 안정적으로 수렴하는 AI-REML 알고리즘을 구현할 수 있게 된다.

Gilmour의 논문에 나오는 observed information matrix, expected information matrix 및 average information matrix는 다음과 같다.

그리고 Larry R. Schaeffer의 2008년도 Animal Models and Estimation of Variance Components에는 다음과 같이 나온다.

observed information matrix와 expected information matrix는 같은 것을 알 수 있으나, average information matrix는 약간 다름을 알 수 있다. 아마도 Schaeffer가 Gilmour의 논문을 인용하며 둘째 식에서 0.5 를 빼먹은 것으로 보인다.

비록 Schaeffer가 0.5를 빼먹었지만 여기서는 Schaeffer의 표기를 기준으로 simplified average를 구해 본다.

셋째 식만 진짜 average 이고 첫째와 둘째는 trace를 qudratic으로 근사하여 average를 간단하게 만든 식이다.

이전 포스트에서 관측 정보 행렬을 유도했다. 이제 이 관측 정보 행렬의 기댓값인 기대 정보 행렬(expected information matrix)을 유도할 차례이다. 여기서 중요한 것은 quadratic form의 기댓값을 구하는 방법이 적용된다는 점이다. quadratic form은 y'Qy 형태 즉 스칼라 값을 만들어 내는 행렬의 곱의 형태이다. 이전 포스트에서도 설명을 했지만 다시 quadratic form의 기댓값을 구해보고, 관측 정보 행렬에 이 공식을 적용하여 기대 정보 행렬을 구해본다.
Gilmour의 Average Information REML An Efficient Algorithm for Variance Parameter Estimation in Linear Mixed Models 논문에는 observed information matrix와 expected information matrix가 다음과 같다.

그리고 Larry R. Schaeffer의 2008년도 Animal Models and Estimation of Variance Components에는 observed information matrix 과 expected information matrix이 다음과 같이 나와 있다.

 
표기법에 약간의 차이가 있을 뿐 같은 식이다. 여기서는 Schaeffer의 observed information matrix를 기댓값을 취하여 expected information matrix를 구해 본다.

REML 방법으로 분산 성분을 추정하는 세 가지 방법이 있다고 하였다. 1) derivative free, 미분하지 않고 그냥 분산 성분 찾아내기 2) 1차 미분하여 0과 같다고 놓고 분산 성분 찾아내기 3) 2차 미분하여 분산 성분 찾아내기.

2차까지 미분하는 이유는 1차 미분한 식에서 식을 0으로 하는 분산 성분 값을 찾기가 힘들고, 추정한 분산 성분을 이용하여 다시 육종가를 구하고 하는데 많은 round를 필요로 한다는 것이다. 원래 likelihood에서 최대 값을 찾기 위하여 1차 미분을 하였다면, 1차 미분 식을 0으로 만드는 값을 빨리 찾기 위하여 2차 미분을 하는 것이다. (자세한 사항은 나중에)

여기서는 1차 미분한 식을 이용하여 2차 미분한 식을 유도해 보고자 한다. 미분 기호식만 쓴다면 2차 미분은 매우 간단하다. 그래서 Johnson and Thomson의 Restricted maximum likelihood estimation of variance components for univariate animal models using sparse matrix techniques and average information 논문에는 다음과 같이 1차 미분식과 2차 미분식이 간단히 나온다.

그리고 Searle의 Variance Component라는 책에도 다음과 같이 나온다.

그러나 계산을 위해서 경우를 나누어 보면 약간 복잡해 진다. 분산을 효과의 분산과 잔차의 분산으로 나누어 본다면 1) 효과의 분산 - 효과의 분산으로 미분 2) 효과의 분산 - 잔차 분산으로 미분 3) 잔차 분산 - 잔차 분산으로 미분하는 경우로 나누어 볼 수 있다. 그래서 Gilmour의 Average Information REML An Efficient Algorithm for Variance Parameter Estimation in Linear Mixed Models 논문에서는 다음과 같이 2차 미분한 식을 다음과 같이 표현하고 있다. 

그리고 Larry Schaeffer의 Animal Model and Estimation of Variance Components 라는 책에서는 다음과 같이 표현하고 있다.

여기서는 2차 미분식을 추정해 보는데 역시나 용어가 어렵다. 1차 미분한 식을 R. A. Fisher가 사용한 용어를 따라 Scoring function 등으로 부른다고 했는데 2차 미분한 식도 이름이 있다. Hessian Matrix. 그렇다면 우리가 구하는 식을 2차 미분한 식이라고 하였는데 Hessian Matrix를 구하는 것인가? 그렇지 않다. 우리가 구하는 식은 Hessian Matrix에 "-" negative를 붙인 식인데 이것을 observed information matrix (관측 정보 행렬???)이라고 한다. 비슷한 이름이 많고 헷갈린다. observed information matrix는 expected informaion matrix를 구하고 observed information matrix와 expected information matrix의 평균인 average information matrix를 구하는 기초가 된다. 여기서 matrix라는 표현을 쓰는데 효과의 분산이 한 개이고, 잔차 분산이 있는 경우 행렬의 (1,1) 원소는 잔차 분산 - 잔차 분산으로 미분한 경우의 값, (1, 2) 원소는 잔차 분산 - 효과의 분산으로 미분한 경우의 값 (2, 2) 원소는 효과의 분산 - 효과의 분산으로 미분한 경우의 값이다. information matrix의 역행렬을 구해서 빠르게 likelihood를 1차 미분한 식을 0으로 만드는 값을 찾는 방법은 나중에 설명하기로 한다. 암튼 여기서는 위 그림의 식을 유도한다.

notation의 약간의 차이가 있지만 암튼 유도를 마쳤다. 

REML의 likelihood function을 이용하여 분산 성분을 추정하는 방법에는 세 가지가 있다고 하였다. 1) 미분하지 않고 2) 1차 미분하고 3) 2차 미분하고.

앞의 포스트는 미분하지 않고 분산 성분을 추정하는 방법에 대한 포스트 였고, 여기서는 미분을 하고 0으로 놓아, likelihood fucntion을 최대로 하는 분산값을 찾는 방법을 설명한다. 미분 자체는 매우 간단하다. determinant의 로그값의 미분이나, 역행렬의 미분에 대해서는 이미 앞선 글에서 설명한 바 있다. 미분이 문제가 아니라 미분한 식을 계산 가능한 식으로 변환하는 것이 문제다. Projection Matrix P가 장악을 한 공식인데 이게 매우 복잡하여 실제로 계산을 하기에 적당하지 않아서, 이것을 계산이 가능한 적당한 식으로 변환하는 것이 문제다. 암튼 여러 참고 문헌과 AI의 도움으로 유도를 하였다. 여러 참고 문헌 중 S. R. Searle의 Variance Component란 책이 많은 도움이 되었다. 대학원생이라면 한 번 보겠지만 지금은 엄두가 나지 않는다.

100년 전에 R. A. Fisher는 log-likelohood function을 1차 미분한 것을 score function이라고 하였다. "추정량이 얼마나 잘 모수를 반영하는가"를 점수(score)처럼 매겨 평가한다는 개념적 비유가 사용된 것이라고 한다. 그래서 지금부터 어디선가 score function이란 말을 들어도 당황하지 말아야 겠다. 

Expectation-Maximization Algorithm

이제 전에도 이용했던 예제를 이용하여 간단히 분산 성분을 추정해 보자.

90개의 자료를 이용하여 구성한 least squares equation이다. 고정효과 1개, 임의 효과 A와 B가 있으며, 고정 효과의 레벨 수(rank of X)는 2(r(X) = 2), 임의 효과의 A의 레벨 수는 3, 임의 효과 B의 레벨수는 4이다. 전체 관측치의 제곱합(total sum of squares)는 356,000이다. 

임의 효과 A의 분산비 초기값은 10, 임의 효과 B의 분산비 초기값은 5로 하여, 해를 구하고(expectation), 그 해를 기초로 likelihood function을 최대화 하는 분산 성분을 추정(maximization)한다. 그리고 구한 분산 성분을 기초로 분산비를 구한 후, expectation - maximization을 반복한다. 수렴 기준(convergence criteria)을 정하고, 수렴할 때까지 반복한다.

# total sum of squares
yty = 356000
yty

# left hand side
lhs = matrix(c(
    50, 0, 5, 15, 30, 5, 10, 20, 15,
    0, 40, 5, 15, 20, 5, 10, 20, 5,
    5, 5, 10, 0, 0, 2, 3, 4, 1,
    15, 15, 0, 30, 0, 5, 7, 11, 7,
    30, 20, 0, 0, 50, 3, 10, 25, 12,
    5, 5, 2, 5, 3, 10, 0, 0, 0,
    10, 10, 3, 7, 10, 0, 20, 0, 0,
    20, 20, 4, 11, 25, 0, 0, 40, 0,
    15, 5, 1, 7, 12, 0, 0, 0, 20
), nrow = 9, byrow = TRUE)
lhs

# right hand side
rhs = matrix(c(3200, 2380, 580, 1860, 3140, 700, 1320, 2400, 1160), nrow = 9)
rhs

# variance ratio
# alpha_a = 40, alpha_b = 10 일 때
alpha_a = 10
alpha_b = 5

lhs2 = lhs

# 임의효과에 분산비 더하기
diag(lhs2)[3:5] = diag(lhs2)[3:5] + alpha_a
diag(lhs2)[6:9] = diag(lhs2)[6:9] + alpha_b

lhs2

# solution 구하기
C = solve(lhs2)
#sol = solve(lhs2, rhs)
sol = C %*% rhs
sol

# rank
N_rx = 88
N_rx

# y'XbZu
ytXbZu =  t(sol) %*% rhs 
ytXbZu


# trCaa
trCaa = sum(diag(C)[3:5])
trCaa

# trCbb
trCbb = sum(diag(C)[6:9])
trCbb

# sigma_e2 추정
sigma_e2 = (yty - ytXbZu ) / N_rx
sigma_e2

# sigma_A2 추정
sigma_A2 = (t(sol[3:5]) %*% sol[3:5] + trCaa * sigma_e2) / 3
sigma_A2

# sigma_B2 추정
sigma_B2 = (t(sol[6:9]) %*% sol[6:9] + trCbb * sigma_e2) / 4
sigma_B2

# new ratio
alpha_a2 = sigma_e2 / sigma_A2
alpha_a2

alpha_b2 = sigma_e2 / sigma_B2
alpha_b2

실행 결과

> # total sum of squares
> yty = 356000
> yty
[1] 356000
> 
> # left hand side
> lhs = matrix(c(
+     50, 0, 5, 15, 30, 5, 10, 20, 15,
+     0, 40, 5, 15, 20, 5, 10, 20, 5,
+     5, 5, 10, 0, 0, 2, 3, 4, 1,
+     15, 15, 0, 30, 0, 5, 7, 11, 7,
+     30, 20, 0, 0, 50, 3, 10, 25, 12,
+     5, 5, 2, 5, 3, 10, 0, 0, 0,
+     10, 10, 3, 7, 10, 0, 20, 0, 0,
+     20, 20, 4, 11, 25, 0, 0, 40, 0,
+     15, 5, 1, 7, 12, 0, 0, 0, 20
+ ), nrow = 9, byrow = TRUE)
> lhs
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
 [1,]   50    0    5   15   30    5   10   20   15
 [2,]    0   40    5   15   20    5   10   20    5
 [3,]    5    5   10    0    0    2    3    4    1
 [4,]   15   15    0   30    0    5    7   11    7
 [5,]   30   20    0    0   50    3   10   25   12
 [6,]    5    5    2    5    3   10    0    0    0
 [7,]   10   10    3    7   10    0   20    0    0
 [8,]   20   20    4   11   25    0    0   40    0
 [9,]   15    5    1    7   12    0    0    0   20
> 
> # right hand side
> rhs = matrix(c(3200, 2380, 580, 1860, 3140, 700, 1320, 2400, 1160), nrow = 9)
> rhs
      [,1]
 [1,] 3200
 [2,] 2380
 [3,]  580
 [4,] 1860
 [5,] 3140
 [6,]  700
 [7,] 1320
 [8,] 2400
 [9,] 1160
> 
> # variance ratio
> # alpha_a = 40, alpha_b = 10 일 때
> alpha_a = 10
> alpha_b = 5
> 
> lhs2 = lhs
> 
> # 임의효과에 분산비 더하기
> diag(lhs2)[3:5] = diag(lhs2)[3:5] + alpha_a
> diag(lhs2)[6:9] = diag(lhs2)[6:9] + alpha_b
> 
> lhs2
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
 [1,]   50    0    5   15   30    5   10   20   15
 [2,]    0   40    5   15   20    5   10   20    5
 [3,]    5    5   20    0    0    2    3    4    1
 [4,]   15   15    0   40    0    5    7   11    7
 [5,]   30   20    0    0   60    3   10   25   12
 [6,]    5    5    2    5    3   15    0    0    0
 [7,]   10   10    3    7   10    0   25    0    0
 [8,]   20   20    4   11   25    0    0   45    0
 [9,]   15    5    1    7   12    0    0    0   25
> 
> # solution 구하기
> C = solve(lhs2)
> #sol = solve(lhs2, rhs)
> sol = C %*% rhs
> sol
            [,1]
 [1,] 64.6313096
 [2,] 59.4225059
 [3,] -2.1363094
 [4,]  0.4994633
 [5,]  1.6368461
 [6,]  5.1063791
 [7,]  2.6402428
 [8,] -2.6432515
 [9,] -5.1033704
> 
> # rank
> N_rx = 88
> N_rx
[1] 88
> 
> # y'XbZu
> ytXbZu =  t(sol) %*% rhs 
> ytXbZu
         [,1]
[1,] 347871.3
> 
> 
> # trCaa
> trCaa = sum(diag(C)[3:5])
> trCaa
[1] 0.1649312
> 
> # trCbb
> trCbb = sum(diag(C)[6:9])
> trCbb
[1] 0.3309886
> 
> # sigma_e2 추정
> sigma_e2 = (yty - ytXbZu ) / N_rx
> sigma_e2
         [,1]
[1,] 92.37198
> 
> # sigma_A2 추정
> sigma_A2 = (t(sol[3:5]) %*% sol[3:5] + trCaa * sigma_e2) / 3
> sigma_A2
         [,1]
[1,] 7.575855
> 
> # sigma_B2 추정
> sigma_B2 = (t(sol[6:9]) %*% sol[6:9] + trCbb * sigma_e2) / 4
> sigma_B2
         [,1]
[1,] 24.16281
> 
> # new ratio
> alpha_a2 = sigma_e2 / sigma_A2
> alpha_a2
         [,1]
[1,] 12.19294
> 
> alpha_b2 = sigma_e2 / sigma_B2
> alpha_b2
         [,1]
[1,] 3.822899

계속 반복하면 임의 효과 A의 분산비는 35.7558로, 임의 효과 B의 분산비는 3.0101로 수렴한다.

위에서는 분산 성분을 추정할 때 추정하려는 효과의 분산 공분산 행렬을 ZZ'로 가정을 했으나 가축 육종의 Animal Model에서는 A로 가정한다. A는 Numerator Relationship Matrix이다. 이것을 한국어로 혈연 계수 행렬로 번역을 하고 있으나 이는 사실과 다르다. 정확히는 혈연 계수의 분자의 행렬이다. A로부터 혈연 계수를 구하려면 off diagonal의 숫자에 두 개체의 근교 계수의 root 값으로 나누어 주어야 한다. 그러나 너무 말이 복잡하므로 그냥 혈연 계수 행렬로 번역을 하는데 정확히 알아두는게 좋을 것 같다. Animal Model에서 개체 효과의 분산 공분산 행렬을 A라고 하였으므로 그것을 반영한 분산 성분 추정이어야 한다. 위에서 잔차의 분산은 변함이 없고 오로지 개체의 효과만 약간의 수정이 있다. 만일 repeatabiliy model 일 경우 영구 환경 효과 P를 고려하는데 MME와 개체의 분산 성분 추정을 위한 식은 다음과 같아 진다.

개체의 분산 성분을 추정할 때 A의 역행렬이 동원되는 것을 알 수 있다. 컴퓨팅 자원이 부족한 시절, A의 역행렬을 구하는 것은 매우 어려운 일이었다. 그런데 C. R. Henderson이 혈통을 보고 바로 A의 역행렬을 구하는 방법을 제시했다. MME를 제시하고, A의 역행렬을 바로 구할 수 있게 됨으로써 육종학의 새로운 지평이 열린 것이다.