Maximum Likelihood Estimation방법으로 혼합 모형의 분산 성분 구하기

자료로부터 Maximum likelihood(ML) 방법으로 평균과 분산을 구해본 적이 있다.

2020.04.24 - [Animal Breeding/Animal Breeding Etc] - Maximum Likelihood Estimation(MLE)

Maximum Likelihood Estimation(MLE)

 

이제 혼합 모형에서 임의 효과의 분산 성분을 구하는 것을 해 보자. 그런데 이전보다 먼저 이해해야 할 것이 참 많다. 다변량 정규 분포 함수나 다변량 우도 함수(likelihood function), 혼합 모형 관측치의 분산-공분산 행렬과 그 역행렬, 공분산 행렬과 MME의 계수 행렬의 행렬식, 투영 행렬(Projection matrix), 그리고 우도 함수의 로그를 취하여 미분하기, 그에 따른 V와 V의 역행렬의 미분, 행렬식의 로그값의 미분, P의 미분, Z'ZTZ'Z의 trace의 변형 등을 먼저 이해하여야 한다. 각각을 간단히 다룬 포스트가 있으므로 미리 공부하면 좋을 것 같다.

이제 혼합 모형에 대한 간단한 설명으로 시작해 보자.

다변량 정규 분포의 likelihood function을 구해본다.

이 함수를 최대화하는 평균과 분산을 구하기 위하여 평균과 분산에 대하여 미분을 하는데 좀 복잡하여 로그를 취한후 미분을 하고 0과 같다고 한다.

등식을 만들고 등식의 좌변과 우변을 간단히 하여 최종적으로 분산 성분을 구할 수 있는 식을 유도한다. 

 

EM(expectation - maximization) Algorithm

육종가(u)와 분산 성분(sigma)은 서로 영향을 미친다. 분산 성분을 알아야 육종가를 구할 수 있고, 육종가를 알아야 분산 성분을 구할 수 있다. 그래서 실제 계산에서는 분산 성분을 모르므로 대충 초기치를 주고, 그 초기치를 기초로 육종가를 구한다(expectation). 그 육종가를 기초로 위 공식에 따라 분산 성분을 구한다(maximization). 이걸 계속 반복해서 안정화 될 때까지 한다. 안정의 기준은 육종가의 변화의 제곱의 합이 얼마 이하(보통 10-12승 에서 10-14승 이하)로 기준을 정하고 계속 반복한다.

그래서 육종가를 (몇 개월 마다 자료가 추가되면) 계속적으로 구해야 하는 경우 분산 성분을 매번 구하면 안 된다. 분산 성분은 고정해 놓고, 자료가 추가된 것에 따라 육종가를 구해야 한다. 그런데 시간이 지나면 육종가도 변하므로 몇 년에 한 번 분산 성분을 재추정하고 다시 몇 년 동안 사용해야 한다. 자료가 많아지고, 개체 수가 많아지고, 동시에 구하는 형질의 수가 증가할 경우(다형질 모형, multiple traits model), 분산 성분 추정을 위하여 컴퓨터가 몇 개월씩 계산할 수도 있음을 알아야 한다.

 

간단한 예제를 통하여 수식을 더 이해해 보자. 아래는 ML 방법으로 임의 효과가 1개인 모형의 분산 성분을 구하는 과정이다. 우선 문제는 다음과 같다.

 

위 least squares equation을 입력하고, 해를 풀고(expectation), 분산 성분을 구해 보자(maximization)

lhs = matrix(c(8, 0, 2, 3, 1, 2, 0,
               0, 14, 3, 3, 2, 4, 2,
               2, 3, 5, 0, 0, 0, 0,
               3, 3, 0, 6, 0, 0, 0,
               1, 2, 0, 0, 3, 0, 0,
               2, 4, 0, 0, 0, 6, 0,
               0, 2, 0, 0, 0, 0, 2
), nrow =7)
lhs

N = sum(diag(lhs)[1:2])
N

rhs = matrix(c(36.3, 67.2, 26.5, 30, 14.1, 26.2, 6.8), nrow = 7)
rhs

Xty = matrix(rhs[1:2], nrow = 2)
Xty

Zty = matrix(rhs[3:7], nrow = 5)
Zty

yty = 500

alpha = 6

lhs2 = lhs
diag(lhs2)[3:7] = diag(lhs2)[3:7] + alpha
lhs2

sol = solve(lhs2, rhs)
sol

b = matrix(sol[1:2], nrow = 2)
b

u = matrix(sol[3:7], nrow = 5)
u

ete = yty - t(b) %*% Xty - t(u) %*% Zty
ete

utu = t(u) %*% u
utu

sigma_e2 = ete / N
sigma_e2

sigma_s2 = (utu + sum(1/diag(lhs2)[3:7]) * sigma_e2 ) / nrow(u)
sigma_s2

alpha2 = sigma_e2 / sigma_s2
alpha2

 

위 R script를 실행한 결과이다.

> lhs = matrix(c(8, 0, 2, 3, 1, 2, 0,
+                0, 14, 3, 3, 2, 4, 2,
+                2, 3, 5, 0, 0, 0, 0,
+                3, 3, 0, 6, 0, 0, 0,
+                1, 2, 0, 0, 3, 0, 0,
+                2, 4, 0, 0, 0, 6, 0,
+                0, 2, 0, 0, 0, 0, 2
+ ), nrow =7)
> lhs
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    8    0    2    3    1    2    0
[2,]    0   14    3    3    2    4    2
[3,]    2    3    5    0    0    0    0
[4,]    3    3    0    6    0    0    0
[5,]    1    2    0    0    3    0    0
[6,]    2    4    0    0    0    6    0
[7,]    0    2    0    0    0    0    2
> 
> N = sum(diag(lhs)[1:2])
> N
[1] 22
> 
> rhs = matrix(c(36.3, 67.2, 26.5, 30, 14.1, 26.2, 6.8), nrow = 7)
> rhs
     [,1]
[1,] 36.3
[2,] 67.2
[3,] 26.5
[4,] 30.0
[5,] 14.1
[6,] 26.2
[7,]  6.8
> 
> Xty = matrix(rhs[1:2], nrow = 2)
> Xty
     [,1]
[1,] 36.3
[2,] 67.2
> 
> Zty = matrix(rhs[3:7], nrow = 5)
> Zty
     [,1]
[1,] 26.5
[2,] 30.0
[3,] 14.1
[4,] 26.2
[5,]  6.8
> 
> yty = 500
> 
> alpha = 6
> 
> lhs2 = lhs
> diag(lhs2)[3:7] = diag(lhs2)[3:7] + alpha
> lhs2
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    8    0    2    3    1    2    0
[2,]    0   14    3    3    2    4    2
[3,]    2    3   11    0    0    0    0
[4,]    3    3    0   12    0    0    0
[5,]    1    2    0    0    9    0    0
[6,]    2    4    0    0    0   12    0
[7,]    0    2    0    0    0    0    8
> 
> sol = solve(lhs2, rhs)
> sol
            [,1]
[1,]  4.42362421
[2,]  4.78319695
[3,]  0.30028734
[4,]  0.19829471
[5,]  0.01222021
[6,] -0.14833635
[7,] -0.34579924
> 
> b = matrix(sol[1:2], nrow = 2)
> b
         [,1]
[1,] 4.423624
[2,] 4.783197
> 
> u = matrix(sol[3:7], nrow = 5)
> u
            [,1]
[1,]  0.30028734
[2,]  0.19829471
[3,]  0.01222021
[4,] -0.14833635
[5,] -0.34579924
> 
> ete = yty - t(b) %*% Xty - t(u) %*% Zty
> ete
         [,1]
[1,] 10.15069
> 
> utu = t(u) %*% u
> utu
          [,1]
[1,] 0.2712234
> 
> sigma_e2 = ete / N
> sigma_e2
          [,1]
[1,] 0.4613951
> 
> sigma_s2 = (utu + sum(1/diag(lhs2)[3:7]) * sigma_e2 ) / nrow(u)
> sigma_s2
           [,1]
[1,] 0.09980162
> 
> alpha2 = sigma_e2 / sigma_s2
> alpha2
         [,1]
[1,] 4.623122
>

 

이건 1 round만 한 것이고, 육종가 또는 분산 성분이 어느 정도 변하지 않을 때까지 계속 반복적으로 실행한다.