REML의 우도함수를 최대화하는 분산을 찾는 방법에는 세 가지가 있다고 이전 글에서 얘기했다. 그 세 방법 중 첫째 방법인 Derivative Free 방법으로 분산을 찾는 방법에 대해서 알아 보자.
Derivative Free라는 것은 미분을 하지 않는다는 것이다. ML의 우도함수를 미분해서 0과 같다고 한 후, 평균과 분산을 구하면 우도 함수를 최대(최소)화하는 것이라 했다. 그런데 미분을 하지 않고 어떻게 찾는다는 것인가? 이 값 저 값을 넣어봐서 최대가 되는지 보고 더 이상 최대가 안 되는 값을 우도 함수를 최대화하는 값이라 하겠다는 것이다. 수학이 꼭 뭘 풀어서 답을 찾는다고 생각할 수도 있지만 너무 복잡해서 풀기 힘들 때는 대충 넣어 보면서 구할 수도 있다. 진짜 최대화하는 값이 아닐 수도 있지만 적당히 비슷하면 세상의 문제를 푸는데 아무 문제가 없기 때문이다. 예를 들어 y = (x - 1)^2 을 최소화하는 x를 찾는다고 했을 때 -10 에서 +10 사이의 값을 대충 넣어 보면 0에서 2 사이에서 최소화되는 것을 알 수 있다. 그러면 범위를 좀 더 좁혀서 0에서 2 사이의 값 10개 정도를 넣어 봐서 범위를 줄이고, 이거를 계속 반복하면 1 또는 1 근처의 값을 구할 수 있다. 물론 이렇게 막무가내로 하지는 않겠지만, 대충 이 값 저 값을 넣어서 원하는 값을 빨리 찾는 법은 많은 수학자들이 수 많은 방법을 연구했다. (대학교때 배운 것도 같지만 나는 잘 모른다.) 암튼 그렇다. 뒤에 예제를 보면 이해가 간다.
우도 함수의 여러 가지 변형 형태에 대해서 얘기했었다. 우도함수의 마지막 항 y'K'(K'VK)-1K'y = y'Py = (y - Xb^)' V-1 ( y - Xb^) 임을 설명하였다. 여기서는 또 다른 변형 형태에 대해서 설명하고 어떻게 적용하는지 알아보자.
이제 다음 예제를 보자.
F는 고정효과이고, A와 B는 임의효과이다. A와 B의 분산을 구해야 한다. A와 B의 분산 값으로 여러 가지 값을 넣어 봐서 어떤 값이 우도함수를 최대화하는지 보면 되는데, 동시에 2개를 처리하려면 경우의 가짓수가 너무 많아지므로 먼저 B의 분산을 고정하고 고정된 B 값에서 어떤 A가 우도함수를 최대화하는지 찾아볼 것이다. 예상했듯이 A가 나오면 구한 A 값으로 고정하고 B의 값을 이 값 저 값을 넣어 봐서 우도함수를 최대화하는 B를 찾는다. 다시 B를 고정하고 A를 찾는다. 이걸 계속 반복해서 안정화될 때까지 한다.
여기서는 분산 값을 찾지 않고 분산비를 찾을 텐데, 잔차 분산을 구하므로 마지막에 각 효과의 분산도 구할 수 있다. 임의 효과 B의 분산비를 10으로 고정하고, 임의 효과 A의 분산비가 5, 10, 20, 30, 40일 때의 우도 값을 구하면 다음과 같다.
여기서는 A의 분산비가 40일 때 우도함수값을 구하는 과정을 설명한다. 그러면 나머지는 알아서 구할 수 있을 것이다.
# 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 일 때, L4 구하기
alpha_a = 40
alpha_b = 10
lhs2 = lhs
# 임의효과에 분산비 더하기
diag(lhs2)[3:5] = diag(lhs2)[3:5] + alpha_a
diag(lhs2)[6:9] = diag(lhs2)[6:9] + alpha_b
lhs2
# solution 구하기
sol = solve(lhs2, rhs)
sol
# rank
N_rx = 88
N_rx
# 잔차 분산 구하기
sigma_e2 = (yty - t(sol) %*% rhs ) / N_rx
sigma_e2
# 잔차 분산 로그값
log_se2 = log(sigma_e2)
log_se2
# y'Py
ytPy = (yty - t(sol) %*% rhs ) / sigma_e2
ytPy
# log of determinant of C*
logC = log(det(lhs2))
logC
# levels of random effect a
q_a = 3
q_a
# no. X log of alpha_a
qlog_a = q_a * log(alpha_a)
qlog_a
# levels of random effect b
q_b = 4
# no. X log of alpha_b
qlog_b = q_b * log(alpha_b)
qlog_b
# 최종 log likelihood
L4 = -0.5 * (N_rx * log_se2 - qlog_a - qlog_b + logC + ytPy)
L4
실행 결과
> # 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 일 때, L4 구하기
> alpha_a = 40
> alpha_b = 10
>
> 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 50 0 0 2 3 4 1
[4,] 15 15 0 70 0 5 7 11 7
[5,] 30 20 0 0 90 3 10 25 12
[6,] 5 5 2 5 3 20 0 0 0
[7,] 10 10 3 7 10 0 30 0 0
[8,] 20 20 4 11 25 0 0 50 0
[9,] 15 5 1 7 12 0 0 0 30
>
> # solution 구하기
> sol = solve(lhs2, rhs)
> sol
[,1]
[1,] 64.8006682
[2,] 59.6949940
[3,] -0.8899764
[4,] 0.1250204
[5,] 0.7649560
[6,] 3.8190836
[7,] 2.3062868
[8,] -2.1370493
[9,] -3.9883211
>
> # rank
> N_rx = 88
> N_rx
[1] 88
>
> # 잔차 분산 구하기
> sigma_e2 = (yty - t(sol) %*% rhs ) / N_rx
> sigma_e2
[,1]
[1,] 96.39973
>
> # 잔차 분산 로그값
> log_se2 = log(sigma_e2)
> log_se2
[,1]
[1,] 4.568503
>
> # y'Py
> ytPy = (yty - t(sol) %*% rhs ) / sigma_e2
> ytPy
[,1]
[1,] 88
>
> # log of determinant of C*
> logC = log(det(lhs2))
> logC
[1] 32.05245
>
> # levels of random effect a
> q_a = 3
> q_a
[1] 3
>
> # no. X log of alpha_a
> qlog_a = q_a * log(alpha_a)
> qlog_a
[1] 11.06664
>
> # levels of random effect b
> q_b = 4
>
> # no. X log of alpha_b
> qlog_b = q_b * log(alpha_b)
> qlog_b
[1] 9.21034
>
> # 최종 log likelihood
> L4 = -0.5 * (N_rx * log_se2 - qlog_a - qlog_b + logC + ytPy)
> L4
[,1]
[1,] -250.9019
이렇게 임의효과 A의 값이 5, 10, 20, 30, 40일 때의 우도함수값을 구하였는데 과연 A값이 얼마일 때 우도함수가 최대가 되는 것일까? 이걸 알려면 2차항 회귀분석을 하면 된다.
# alpha_b = 10으로 고정하고
# alpha_a가 5, 10, 20, 30, 40일 때의 L4
# L4를 최대화하는 alpha_a 구하기
# 2차 회귀 곡선으로 적합
Q = matrix(c(
1, 5, 25,
1, 10, 100,
1, 20, 400,
1, 30, 900,
1, 40, 1600
), nrow = 5, byrow = TRUE)
Q
Y = matrix(c(
-251.4442,
-251.1504,
-250.9822,
-250.9274,
-250.9019
), nrow = 5)
Y
beta = solve(t(Q) %*% Q, t(Q) %*% Y)
beta
m_alpha_a = - (beta[2] / (2 * beta[3]))
m_alpha_a실행 결과
> # alpha_b = 10으로 고정하고
> # alpha_a가 5, 10, 20, 30, 40일 때의 L4
> # L4를 최대화하는 alpha_a 구하기
> # 2차 회귀 곡선으로 적합
> Q = matrix(c(
+ 1, 5, 25,
+ 1, 10, 100,
+ 1, 20, 400,
+ 1, 30, 900,
+ 1, 40, 1600
+ ), nrow = 5, byrow = TRUE)
> Q
[,1] [,2] [,3]
[1,] 1 5 25
[2,] 1 10 100
[3,] 1 20 400
[4,] 1 30 900
[5,] 1 40 1600
>
> Y = matrix(c(
+ -251.4442,
+ -251.1504,
+ -250.9822,
+ -250.9274,
+ -250.9019
+ ), nrow = 5)
> Y
[,1]
[1,] -251.4442
[2,] -251.1504
[3,] -250.9822
[4,] -250.9274
[5,] -250.9019
>
> beta = solve(t(Q) %*% Q, t(Q) %*% Y)
> beta
[,1]
[1,] -2.516015e+02
[2,] 4.488460e-02
[3,] -6.979310e-04
>
> m_alpha_a = - (beta[2] / (2 * beta[3]))
> m_alpha_a
[1] 32.15547
즉 B의 분산비가 10일 때 우도함수를 간단히 다음과 같은 2차 함수로 나타낼 수 있다는 뜻이다.
y = -0.000697931 x^2 + 0.0448846 x - 251.6015, 여기서 x는 임의효과 A의 분산비
이 2차 함수를 최대화하는 x는 - B/ 2A 이므로 0.0448846 / (2 * 0.000697931) = 32.15547 이 된다. 이제 임의효과 A의 분산비를 32.15547로 고정하고, 임의 효과 B의 값을 2, 3, 4, ..., 10로 바꾸어 가며 우도함수 값을 구하고 위와 같이 2차항 회귀분석을 하여 임의 효과 B의 값을 정한다. 다시 A의 값을 32 전후로 바꾸어 가며 A 값을 바꾼다. 이때 점점 값 사이의 간격이 작아질 것이고 결국 어떤 수준 이하로 수렴할 것이다.
그런데 여기서 우리가 구한 값이 global maximum이란 보장이 없다. 여러 개의 산이 있는데 비교적 높이가 낮은 산 근처에서 시작을 해서 눈 감고 찾다 보니 높이가 낮은 산의 높이를 구할 수도 있다. 이걸 local maximum이라 한다. 암튼 전체 구역에서 최고의 높이를 보이는 산을 찾으려면 여러번 여러 군데에서 시작을 해야 한다. 그래도 결과가 같게 나오면 그걸 global maximum이라 인정한다.