저작자표시 - 비영리 - 변경금지 2.0 대한민국 이용자는아래의조건을따르는경우에한하여자유롭게 이저작물을복제, 배포, 전송, 전시, 공연및방송할수있습니다. 다음과같은조건을따라야합니다 : 저작자표시. 귀하는원저작자를표시하여야합니다. 비영리. 귀하는이저작물을영리목적으로이용할

저작자표시 - 비영리 - 변경금지 2.0 대한민국 이용자는아래의조건을따르는경우에한하여자유롭게 이저작물을복제, 배포, 전송, 전시, 공연및방송할수있습니다. 다음과같은조건을따라야합니다 : 저작자표시. 귀하는원저작자를표시하여야합니다. 비영리. 귀하는이저작물을영리목적으로이용할수없습니다. 변경금지. 귀하는이저작물을개작, 변형또는가공할수없습니다. 귀하는, 이저작물의재이용이나배포의경우, 이저작물에적용된이용허락조건을명확하게나타내어야합니다. 저작권자로부터별도의허가를받으면이러한조건들은적용되지않습니다. 저작권법에따른이용자의권리는위의내용에의하여영향을받지않습니다. 이것은이용허락규약 (Legal Code) 을이해하기쉽게요약한것입니다. Disclaimer

이학석사학위논문 디리크레과정혼합모델을이용한유전자자료분석 THE APPLICATION OF DIRICHLET PROCESS MIXTURE MODEL TO DNA MICROARRAY DATA 2017 년 8 월 서울대학교대학원 통계학과 박정연

국문초록 군집화 분석에서는 군집의 개수가 정해져 있지 않는 경우가 대부분이다. 이 러한 경우에는 보통 디리크레 과정을 이용하여 추론한다. 혼합모형에 디리크레 과정을 적용시켜 다양한 군집화 분석 모형을 연출 할 수 있다. 본 논문에서는 우선 유전자 자료의 군집화 분석을 위한 이론 배경에 대하여 소개하고, 디리크레 과정과 무작위 효과 혼합 모형을 결합시켜 유전자 자료의 군집 분석에 적용 해 본다. 이론 부분에서는 디리크레 과정과 디리크레 과정 무작위 효과 혼합 모형, 그리고 추론에서 사용되는 표본 추출법 메트로폴리스-헤이스팅스 알고리즘을 살펴 보도록 한다. 주요어 : 디리클레 과정, 혼합모형, 군집 분석, 유전자 자료 학번 : 2015-22123 i

Contents 1 서론 1 2 이론소개 3 디리크레 과정............................ 3 2.1.1 폴랴 항아리 열....................... 5 2.1.2 세수라만 표현........................ 7 2.2 디리크레 과정 무작위 효과 혼합 모형............... 8 2.3 알고리즘............................... 11 2.1 2.3.1 메트로폴리스-헤이스팅스 알고리즘............ 12 2.3.2 모수 벡터 ξ의 표본추출................... 13 2.3.3 위치 확률 행렬 P 의 추정.................. 15 3 시뮬레이션 자료분석 17 4 유전자 자료분석 21 5 결론 24 참고문헌 25 ii

List of Tables 2.1 군집 Zi 의 표본 추출 법에서 사용 되는 헤이스팅스 비율 값의 계산법 대응 표........................... 13 3.1 p λφ, λτ, λ, 각 군집의 평균과 분산............... 17 iii

List of Figures 3.1 시뮬레이션자료........................... 18 3.2 시뮬레이션자료군집분석결과 ( 평균벡터 )............ 19 3.3 시뮬레이션자료군집분석결과 ( λ φ, λ τ, λ ɛ )......... 20 4.1 유전자자료군집분석결과 ( 평균벡터 ).............. 21 4.2 유전자자료군집분석결과 ( λ φ, λ τ, λ ɛ )............ 22 4.3 유전자자료군집분석결과 (PCA)................. 23 iv

Chapter 1 서론 실제 군집 분석 과정에서 군집의 개수를 예측하는 것은 쉽지 않는 작업이다. 그래서 많은 학자들이 군집의 개수를 미리 지정하지 않고 군집 분석을 할 수 있는 방법을 찾게 되었다. 그 결과로 디리크레 과정, 인디언 뷔페 과정을 이용한 비모수 베이지안 방법이 각광을 받았다. 본 논문에서는 디리크레 과정(Dirichlet process)을 소개한다. 디리크레 과정은 다양한 분포, 모형과 혼합하여 군집 분석 에 널리 사용 되고 있다. 또한, 디리크레 프로세스 혼합 모형에 대하여, 군집의 개수를 추론하는 과정에서 사용되는 조각 추출 법[Walker(2011)], 무공백 추출 법[Neal (2000)], 등 다양한 표본 추출 법이 연구 되었다. 본 논문에서는 붕괴 깁스 표본 추출법(collapsed Gibbs sampling) 중 MCMC(Markov Chain Monte Calo) 메트로폴리스-헤이스팅스(Metropolis Hastings) 알고리즘[Neal (2000)]을 발전시킨 Fu (2013)에서 다룬 새로운 메트로폴리스-헤이스팅스 알고리즘을 살 펴 보았다. 그리고 자료의 군집 분석을 하기 위한 무작위 효과 혼합 모형[Fu (2013)]에 대해서도 소개 하였다. DNA 마이크로어레이 데이터의 분류는 질병의 잠재적인 타입을 파악하고 정 확한 치료 법을 지정하는 데에 유용하다. 본 논문에서 사용되는 자료는 Housden (2011)에서 소개한 초파리 성체 근육세포에 관한 마이크로 어레이 유전자 데이 1

터이다. 초파리의 성체에 5분간의 자극을 주어 활성화 시킨 후 150분 동안 18 개의 불균등한 시점에서 4개 생물학 복제물에 대해 근육세포 mrna 농도를 측정하는 과정을 통해 얻은 163개 차별적 유전자 데이터에 디리크레 프로세스 무작위 효과 혼합 모형을 적용시켜 유전자 표현 패턴을 분류하는 군집 분석을 시행하였다. 본 논문은 총 4 부분으로 나누어졌다. 2 절에서는 디리크레 과정, 디리크레 과 정 무작위 효과 혼합 모형, 그리고 메트로 폴리스-헤이스팅스 알고리즘을 포함한 이론적인 배경 지식을 소개하고, 3절에서는 시뮤레이션 자료에 대한 분석을 다 루며, 4절에서는 유전자 자료에 적용시킨 사례를 살펴 보았다. 마지막으로 이에 대한 결론을 내리며 논문을 마무리 한다. 2

Chapter 2 이론소개 2장에서는 기본 이론 배경에 대해 살펴 보도록 한다. 우선 디리크레 과정 혼합 모형을 소개하기 위해 디리크레 과정에 대해 전반적으로 리뷰하였다(2.1). 2.2에서는 분석에 사용한 기본 모형 무작위 효과 혼합 모형에 대하여 살펴보고, 디리크레 과정 혼합 모형을 어떻게 적용 하였는지를 설명하였다. 마지막으로, 분석과정에서 사용한 MCMC 깁스 표본 추출법 메트로폴리스-헤이스팅스 알고 리즘을 알아 보도록 한다(2.3). 2.1 디리크레 과정 Ferguson, Thomas(1973) 에서 소개된 디리크레 과정은 데이터 마이닝, 기계 학습(machine learning), 등 다양한 분야에서 사용 되고 있다. 디리크레 과정은 일종의 확률 과정으로서 비모수 베이지안 모형에, 특히 디리크레 과정 혼합 모형 에서 많이 사용 되고 있다. 디리크레 과정은 확률 분포에 대해 정의 된 분포이다. 즉, 확률 변수가 어떤 종류의 수가 아닌 확률 분포로 정의 되어 있다는 것이다. 또한, 디리크레 과정에서 추출된 분포들은 이산형이지만 유한 개의 모수로 설명 3

할 수 없다. 그래서, 디리크레 과정을 이용한 군집 분석의 모형은 비모수 모형 이다. 또한, 사후 분포를 추론 할 때 편의상 모수 들의 사전 분포를 폴랴 항아리 열(Polya Urn sequence ) 혹은 세수라만 표현(Sethuraman s representation) 방 식을 사용하기도 한다. 우선 디리크레 과정 부터 살펴 보자. 측도 가능한 공간(measurable space) Θ 상에서 정의 된 확률 측도 G0 를 기본 분포(base distribution)라 하고, 양의 실수 α를 수렴 모수(concentration parameter)라 하자. 여기서 αg0 를 기본 측도(base measure) 라고 한다. Θ의 임 의의 유한한 측도 가능 분할 A1, A2,, Ak 에 대해서, 다음과 같은 조건을 만족 하면, 확률 분포 G가 디리크레 과정 DP를 따른 다(G DP (G0, α))고 할 수 있다. (G(A1 ),, G(Ak )) Dir(αG0 (A1 ),, αg0 (Ak )) (2.1) 확률 분포 G가 무작위적(random)이므로, Θ의 임의의 유한한 측도가능 분할 A1, A2,, Ak 에 대해서, (G(A1 ),, G(Ak ))도 무작위적인 것을 알 수 있다. 또한, DP를 따르는 임의의 확률 분포 G의 주변 확률 분포(marginal distribution) 는 디리크레 분포를 따른 다는 것을 위 식을 통하여 알 수 있다. DP의 존재성은 콜모고로프 일치성 조건(Kolmogorov s Consistency Theorem)을 이용하여 증명 할 수 있다. 확률의 성질과 디리크레 분포의 성질을 이용 하여 다음과 같은 디리 크레 과정에 관한 성질 들을 유도 할 수 있다. 임의의 측도 가능한 집합 A Θ 에 대하여, (1) G(A) Beta(αG0 (A), αg0 (Ac )), (2) DP의 평균은 E[G(A)] = G0 (A)이고, (3) 분산은 V ar[g(a)] = G0 (A)(1 G0 (A))/(α + 1)이며, (4) 만약 무작위적 변수 X R에 대해서, X P P 이면, X G0 이다. 여기서, (1)은 디리크레 분포의 감마표현을 통해 증명 할 수 있고, (2)와 (3)은 평균과 분산의 공식을 이용하여 증명 할 수 있으며, (4)에 대한 증명은 디리크레 과정의 정의와 (1), (2)를 이용 하여 얻을 수 있다. 이에 관한 자세한 증명은 본 4

논문에서 생략 하도록 한다. 다음은 디리크레 과정의 사후 분포에 대해 살펴 본다. 만약 G DP (α, G0 )이고, θ1, θ2,, θn G이면, G의 사후 분포는 여전히 DP이다: N G θ1,, θn DP (α + N, X 1 (αg0 (θ) + δ(θ = θi ))) α+n i=1 (2.2) 여기서, δ(θ = θi ) 는 θ = θi 일 때 1의 값을 갖고, 다른 경우에는 0의 값을 갖는 함수를 뜻한다. 측도 가능한 공간 Θ 상의 임의의 분포에 대하여 DP가 결합 사전 분포(conjugate prior)임을 알 수 있다. 사후 분포를 추정하는 과정에서 계산의 편의를 위하여 보통 사전 분포의 디리 크레 과정을 변형하여 사용한다. 2.1, 2.2절에서 흔히 사용되는 폴랴 항아리 열, 그리고 세수라만 표현(막대 자르기 표현)을 살펴 보도록 한다. 2.1.1 폴랴 항아리 열 폴랴 항아리 열을 소개하기 위해 우선 다음과 같은 상황을 생각해 보자. 가령 항아리에서 k개의 색깔이 있고, 이를 i = 1,, k와 같이 표시 하면, i 색을 가진 P 공의 개수는 αi 이고, A = ki=1 αi 이다. 이 항아리에서 하나의 공을 뽑아서, 그 공의 색을 xj 라고 하자. 또한, xj 는 j 번째 뽑은 공의 색이며, j = 1, 2, 이다. 규칙은 다음과 같다. 항아리에서 공 하나를 뽑고, 그와 같은 색의 공을 뽑은 공과 같이 다시 항아리에 넣는다. 즉, 하나의 빨간 공이 뽑힌 경우에, 두 개의 빨간 색의 공을 다시 항아리에 넣는 뜻이다. 이런 규칙에 따라서 j = 1, 2,, 번째 뽑은 공이 i (i = 1,, k)색을 가진 확률을 계산해 보면 다음과 같다. P αi + nj=1 δxj (i) P (Xn = i X1,, Xn 1 ) = A+n 5 (2.3)

여기서, δxj (i)는 j번째 뽑은 공의 색이 i이면 1이고, 아니면 0이다. 식 (2.3)을 통해서 Xj 들의 결합분포를 다음과 같이 구할 수 있다. Qk P (X1 = x1,, Xn = xn ) = αi (αi + 1) (αi + ni 1) A(A + 1) (A + n 1) i=1 Qk = 여기서, ni = Pn j=1 i=1 [αi ]ni (2.4) [A]n I(Xj = i)는 i색의 공을 뽑은 휫수를 뜻한다. 또한, 이 확 률은 xj 들의 순서에 영향을 받지 않으므로, X1, X2, 는 서로 교환 가능하다. 디피네티 정리를 적용시키면: i.i.d. P = (P1,, Pk ), s.t.x1, X2, P P (2.5) P 의 분포를 알아보기 위하여 우선 적률에 대해서 계산해보았다. P n = ki=1 ni 이라고 하면, E(pn1 1 pnk k ) = E[P rob(x1 = 1,, Xn1 = 1, Xn1 +1 = 2,, Xn1 +n2 = 2,, Xn1 + +nk = k P )] = P rob(x1 = 1,, Xn1 = 1, Xn1 +1 = 2,, Xn1 +n2 = 2,, Xn1 + +nk = k) Qk [αi ]ni = i=1 (2.6) [A]n 또한, 감마 표현을 이용하여, 디리크레 과정을 따르는 P = (P1,, Pk ) 의 적 룰을 구해 보면, E[pn1 1 pnk k ] Qk = i=1 [αi ]ni [A]n (2.7) 이다. 식 (2.6)과 식(2.7)이 동일하다는 것을 알 수 있다. 그러므로 폴랴 항아리 열의 디피네티 측도는 디리크레 과정임을 알 수 있다. 본 논문에서는 대표적으로 이산형인 경우만 다루었지만, 연속형인 경우도 유사한 과정을 통해 같은 결론을 얻을 수 있다. 6

폴랴 항아리 열과 같이 자주 등장 하는 디리크레 과정 사전 분포의 또 다른 표현 방식으로 세수라만 표현이 있다. 세수라만 표현에 관한 내용은 다음 절에서 다루었다. 2.1.2 세수라만 표현 세수라만 표현에 접근 하기 위하여 세수라만 구조(Sethuraman s construction) - 막대 자르기 구조(Stick - Breaking construction)라고도 부르다 - 을 통 해서 살펴 보도록 한다. 길이가 1인 막대를 연속으로 k번 자르는 상황을 생각해 보자. k = 1일 때, 임의의 비율 β1 < 1를 정하여 막대를 자르면, 길이가 π1 = β1 의 막대를 얻을 수 있다. k = 2일 때, 임의의 비율 β2 를 정하여, 남은 1 β1 의 막대를 자르면, 길이가 π2 = β2 (1 β1 )의 막대를 얻을 수 있다. k가 커질 수록 남은 막대의 길이는 더 짧아지는 것을 쉽게 알 수 있다. 이 과정을 식으로 나타나면 다음과 같다: βk Beta(1, α) πk = βk (2.8) k 1 Y (1 βl ) (2.9) l=1 남은 막대의 길이의 집합은 수열 {πk } k=1 로 표현 할 수 있다. 남은 막대의 길이의 분포는 수렴 모수 α에 의하여 결정 된다. α가 크면, 분포는 균일 분포(Uniform distribution)에 근사하다. 베타분포의 성질에 의하여 E(βk ) = 1 임을 1+α 알 수 있다. 이를 통해서, α가 작은 값을 가질 때, βk 는 큰 값을 취하게 되고, 남은 막대의 길이는 급속히 감소하다는 것을 알 수 있다. ηk G0 일 때, 무작위 이산 확률 분포 G(θ) = X πk δ(θ = ηk ) k=1 7 (2.10)

는 DP (α, G0 )를 따르다는 것을 알 수 있다.[증명과정은 A Simple Proof of the Stick-Breaking Construction of the Dirichlet Process, John Paisley, Department of Computer Science Princeton University, Princeton, NJ, jpaisley@princeton.edu 참조] 디리크레 과정에서 무작위하게 추출된 분포는 모두 식 (2.10)과 같은 식으로 표현 가능하다. 이러한 표현 방식을 세수라만 표현이 라고 한다. 또한, 식(2.9) -세수라만 표현- 에 의하여, 디리크레 과정에서 추출한 분포 들은 무한 이산 분포(infinite discrete distribution)라는 것을 알 수 있다. 앞서 2.1절에서 디리크레 과정에 대하여 간략하게 소개를 하였다. 2.2절에서는 디리크레 과정을 이용한 디리크레 과정 무작위 효과 혼합 모형에 대하여 살펴 보도록 한다. 2.2, 2.3절에서 다루는 내용은 [Fu (2013)]에서 다룬 내용에 대한 간략한 리뷰 소개이다. 2.2 디리크레 과정 무작위 효과 혼합 모형 우선 무작위 효과 혼합 모형에 대해 소개를 하고, 모수 예측하는 과정에서 사용 되는 디리크레 과정 혼합 모형에 대해 살펴 보도록 한다. J개의 시점(시점들 사이의 간격은 규칙적이지 않다.)에서, N개의 객체(object) 를 R번 반복하여 측정을 했다고 가정하면, Mijr, i = 1,, N, j = 1,, J, r = 1,, R을 j번째 시점에서, i번째 객체의 r번째 반복 측정의 관측 값이라고 하자. 데이터를 다음과 같이 벡터의 형식 표현하고, K개의 혼합 요소(mixture components) - 보통 여러 개의 혼합 요소가 하나의 군집(cluster)을 정의 할 수 있지만, 여기에서는 혼합 요소와 군집을 동등한 개념으로 취급 하도록 한다 - 가 있다고 가정을 하면, Mi = (Mi11,, Mi1R,, MiJ1,, MiJR )T (2.11) 이고, Mi 들은 서로 독립하고, K개의 혼합 요소를 가진 혼합 비율(mixing proportion)이 wk 인 혼합 분포를 따른다. 만약, gk 를 k번째 혼합 요소의 확률 밀도 8

함수, 평균 벡터 Θ = (Θ1,, ΘK )와 공분산 행렬 Σ = (Σ1,, ΣK )를 혼합 분포의 모수라고 하면, Mi 의 확률 밀도 함수는, f (Mi Θ, σ) = K X wk gk (Mi Θk, Σk ) (2.12) k=1 이다. w를 혼합 비율들의 집합, Zi 를 i번째 객체가 속하는 군집이라고 하면, P rob(zi = k w, Θ, Σ) = wk (2.13) 이다. 또한 자료 M = (M1,, MN )가 주어 졌을 때, 군집의 사후 확률을 P rob(zi = k M ), i = 1,, N, k = 1,, K라 하고, 이를 행렬 P (N K) 로 표현한다. 행렬 P 를 추정하는 것이 목적이며, 추정된 P 를 통해, 군집 분석 을 하고자 한다. 만약, i 번째 객체가 k 번째 군집에서 추출 되었다고 했을 때, Zi = k, 무작위 효과 모형은, Mijr {Zi = k} = Θk + φki + τijk + kijr, (2.14) E(Mijr {Zi = k}) = Θkj, i.i.d. φki {Zi = k, λkφ } N (0, λkφ ), i.i.d. τijk {Zi = k, λkτ } N (0, λkτ ), (2.15) i.i.d. kijr {Zi = k, λk } N (0, λk ). 이다. 여기서, Θkj : j시점에서의 고정 효과(참값), φki : 군집 안의 무작위 효과, τijk : 교호작용 무작위 효과, kijr : 독립적 효과, 이다. i번째 객체가 속하는 군집이 주어졌을 때, Mi 는, i.i.d. Mi {Zi = k, Θk, λkφ, λkτ, λkτ } NJR (Θkagg, Σkagg ), 9 (2.16)

T T 즉, 평균 벡터 Θkagg = (Θk,, Θk )T (Θk 를 R번 반복), 공분산 행렬 Σkagg = [Cov(Mijr, Mij 0 r0 )]JR JR 인 분포를 따른다. 공분산 행렬의 각 요소를 구하는 식 은 다음과 같다. Cov(Mijr, Mij 0 r0 ) = λkφ + λkτ I(j = j 0 ) + λk I(j = j, r = r0 ). (2.17) 이 모형에서 우리가 관심 있는 요소 들은 군집의 개수 K, 구성 모수 들 Θk, λkφ, λkτ, λk, 그리고 사후 할당 확률(posterior allocation probability) 행렬 P (N K)이다. 관심 있는 요소 들에 관한 추정을 하기 위해서 이 모형에 디리크레 과정을 적용 하여 분석에 사용 될 모형을 도출 할 것이다. 다음 부분에서 이에 관한 내용을 살펴 보도록 한다. 앞서 소개한 바와 같이 디리크레 과정은 군집의 개수 K에 대한 사전 정보 없이 관측치와 모수 공간에서의 추출(sampling)에 의존하여 군집화하는 방법이다. 만 약 F 가 다변량 정규 분포(multivariate normal distribution)이고, Mi F (γi ), γi 는 모수들의 집합, γi = {Θi, λφ, λτ, λ } 이라하고, γi 는 수렴 모수가 α, 기본 측 도가 G0 인 디리크레 과정을 따르는 무작위 분포 G를 따른다고 가정하면, 모형을 다음과 같이 표현 할 수 있다: γi G (2.18) G DP (α, G0 ), α 0. (2.19) 위 식에서 γi 들은 같은 분포를 따르지만, 서로 독립이 아님을 알 수 있다. 다음은 모수 γi 들을 디르크레 과정을 이용하여 군집분석하는 과정에 대하여 살펴 보도 록 한다. 군집분석의 과정을 설명하기 위해 우선 G0 에서 γ1 을 추출한다. 그다음 G0 에서 추출한 γ2 의 값은 두 가지 경우가 있다. 1/(1 + α)의 확률로 γ1 과 같은 값을 갖는 경우와 α/(1 + α)의 확률로 γ1 과 다른 값을 갖는 경우이다. 이 과정을 일반화 시켜서 n개의 값을 추출하고, n + 1번째 추출된 값은 다음과 같은 분포를 따른다 10

[Antoniak(1974)]: P r(γn+1 = γ γ1,, γn, α) = P n i=1 1(γi =γ), γ {γ1,, γn }, n+α α, n+α (2.20) γ / {γ1,, γn }. 이 분포에서는 무작위적 군집 개수 K를 포함한 N 개의 값 γ1,, γn 의 분할에 대한 분포 정보가 포함 되어 있다. 수렴 모수 α가 주어졌을 때, 각 γi 들과 대응 되는 군집 Zi, i = 1,, N 과 군집의 개수 K의 분포는 다음과 같다 (Nl : l번째 군집의 크기): K Y Γ(α) αk P r(z1,, ZN, K α > 0) = (Nl 1)!, Γ(α + N ) l=1 P r(z1 = = ZN, K = 1 α = 0) = 1. (2.21) (2.22) 본 논문에서는 위의 분포를 사전 분포(DP (α))로 사용하기로 한다. 식 (2.20)은 폴랴 항아리 열 표현을 이용한 중국집 과정이며, Neal(2000)에서 다룬 메트로폴 리스-헤이스팅스 표본 추출법에서 사용된 표현이다. 여기서 수렴 모수 α가 큰 값을 가졌을 때, γi 들은 높은 확률로 새로운 값을 갖게 되고; 반대로 수렴 모수 α의 값이 작을 수록, γi 들은 원래 존재 했던 값을 갖는 확률이 더 높다. 때문에 α의 값을 이용하여 군집의 폭을 효과적으로 통제할 수 있다. 또한, 식 (2.20)을 이용하여 각 군집에 분류된 원소들의 정보를 얻을 수 있다는 것을 알 수 있다. 2.3에서는 Neal(2000)에서 소개된 내용을 바탕으로 이번 모형에 적합한 진화된 메트로폴리스-헤이스팅 표본 추출법에 대해서 다룰 것이다. 2.3 알고리즘 2.2에서 소개한 디리크레 과정 사전 분포를 따르는 군집 Zi 의 분할 과정 은 MCMC 알고리즘의 핵심이다. 2.3에서는 비결합(nonconjugate) 사전 분포를 따르는 모수와 교차 효과가 있는 모형에도 적합한 메트로폴리스-헤이스팅스 알 고리즘에 대해 소개 할 것이다. 11

2.3.1 메트로폴리스-헤이스팅스 알고리즘 메트로폴리스-헤이스팅스 알고리즘은 MCMC 알고리즘에서 군집 Zi 의 업데 이트 과정에서 사용된다. Zi 의 현재의 값이 z 0 이고, 군집의 개수 K는 k 0 라고 가 정한다. 다음 단계의 업데이트하는 과정에서는 Zi 에게 새로운 값 z 를 추천하고, 이에 따른 군집의 개수는 k 이다. 또한, 디리크레 과정 무작위 효과 모형에서의 관심 모수들의 집합을, 1 K 1 K ξ = {K, Θ1,, ΘK, λ1φ,, λk φ, λτ,, λτ, λ,, λ, Z1,, ZN, α} 로 표현한다. 추천한 새로운 값을 선택할 확률은 min(1, H)이다. ξ에서 Zi 를 제 외한 모든 모수의 현재 추정값을 로, i번째 객체를 제외한 모든 객체의 군집을 z i 로 표현하면, 확률H는 다음과 같이 계산 할 수 있다. π(zi = z ) g(zi = z 0 Zi = z ) H= π(zi = z 0 ) g(zi = z Zi = z 0 ) P r(mi Zi = z, )P r(z1,, z,, zn, k α) g(z 0 z ) = P r(mi Zi = z 0, )P r(z1,, z 0,, zn, k 0 α) g(z z 0 ) P r(mi z, ) P r(z, k z i, α) g(z 0 z ). = P r(mi z 0, ) P r(z 0, k 0 z i, α) g(z z 0 ) (2.23) 여기에서, i번째 객체를 제외한 객체들의 군집 z i 는 i번째 군집을 업데이트 할 때 변하지 않는다. 디리크레 과정 사전 분포를 가정 할 때, 조건부 확률 P r(z 0, k 0 z i, α)를 다음과 같이 계산 할 수 있다(Nz : z번째 군집의 크기): P r(zi = z, K = k Z i = z i, α) Nz 1, Zi 가 단일원소 군집이 아닌 경우, N 1+α = α, Zi 가 단일원소 군집인 경우. n 1+α (2.24) Neal(2000)에서 제시한 MH 알고리즘에서는 제안(proposal) 분포 g가 식(2.24) 와 같은 형태를 가지므로서 보다 간단한 헤이스팅 비율 H = P r(mi z, ) 를 P r(mi z 0, ) 사 용하였다. Fu (2013)에서 다룬 MH 알고리즘은 Neal(2000)에서 다룬 알고리즘 의 객체의 군집 간의 이동이 느린 점을 보완하여 발전시킨 알고리즘이다. Fu 12

(2013)에서 제시한 군집의 제안 분포는 현제 객체가 속한 군집을 제외한 정수열 {1,, k 0 + 1}에서 정의된 이산 균등분포이다. 현제 객체가 속한 군집을 제안 범위에서 제외 시켰기 때문에 객체를 새로운 군집 혹은 크기가 작은 군집으로 더 쉽게 이동 할 수 있다. 그리고 새로운 군집의 채택 여부는 헤이스팅스 비율의 값에 따라 결정이된다. 헤이스팅스 비율은 다음과 같이 계산할 수 있다. Table 2.1: 군집 Zi 의 표본 추출 법에서 사용 되는 헤이스팅스 비율 값의 계산법 대응 표 현재 단일 군집 소속 기존 군집 제안 k k0 1 Yes Yes -1 2 Yes No 0 3 No Yes 0 4 No No 1 헤이스팅스 비율 P r(mi z, ) Nz k0 P r(mi z 0, ) α k0 1 P r(mi z, ) P r(mi z 0, ) P r(mi z, ) Nz P r(mi z 0, ) Nz0 1 P r(mi z, ) α k0 P r(mi z 0, ) Nz0 1 k0 +1 여기서 k 0 와 k 는 각각 i번쨰 객체가 분류된 군집의 변화 전과 후의 전체 군집 개수를 표기한다. i번째 객체가 현제 분류된 군집 Zi = z 0 라 하고, {1,..., z 0 1, z 0 + 1,, k 0 + 1}상에서 정의된 이산 균등분포에서 추출한 제안 군집이 z 라고 가정 하면, 헤이스팅스 비율의 값은 표1에 의하여 결정 된다(supplemental material [Fu(2013)]). 2.3.2 모수 벡터 ξ의 표본추출 모수들의 집합 ξ의 표본추출 과정은 MCMC 알고리즘을 사용하며, 아래와 같은 2 단계를 반복하여 실행한다. 1. Zi (i = 1,, N )를 MH 알고리즘을 이용하여 업데이트 하고, 2. 단계1에서 업데이트한 군집 정보를 조건으로 사용하여 ξ에 속한 다른 모수들 13

을 깁스 또는 MH 알고리즘을 이용하여 업데이트한다. ξ에 속한 다른 모수들의 업데이트 과정을 살펴보기 위해 모형 (2.16)을 계층모 형으로 다시 표현한다. 첫 번째 계층은 j시점에서의 i번째 객체의 반복 측정 관측치에 대한 분포이 다.(µij : j시점에서 i번째 객체의 평균) i.i.d. Mijr {Zi = k, µij, λk } N (µij, λk ), 두 번째 계층은 k번째 군집에 분류된 모든 객체의 평균 벡터에 대한 분포이다.(S: 모든 원소가 1인 J J행렬, I: 대각선에 있는 원소가 1인 J J 단위 행렬.) i.i.d. µi {Zi = k, Θk, λkφ, λkτ, λk } NJ (Θk, λkφ S + λkτ I). 변수 µi 와 모수 Θk, λkφ, λkτ, λk 에 대한 업데이트 알고리즘은 다음과 같다. 1. 현재 k번째 군집에 분류된 각 객체의 평균 벡터 µi (1 J),i = 1,, Nk 에 대한 업데이트를 한다. 평균 벡터는 다변량 정규분포를 따르며, 사후 분포는 다음과 같다. π(µi ) = π(µi Mi, Θk, λkφ, λkτ, λk ) R Y p(mi1r,, MiJr µi, λk I)p(µi Θk, λkφ S + λkτ I) r=1 NJ (µm, ΣM ) 분포 함수를 통해 계산한 결과: ΣM = {R(λk ) 1 I + (λkφ S + λkτ I) 1 } 1, µm = ΣM {(λkφ S + λkτ I) 1 Θk + R(λk ) 1 y m }, R 1X y m = Mi r. R r=1 14

2. 군집 k에 분류된 객체들의 평균 벡터를 이용하여 군집 k의 평균 Θk 에 대한 업데이트를 한다. 3. 군집 k에 분류된 객체들의 모든 데이터를 사용하여 MH 알고리즘을 통해 λkφ, λkτ, λk 를 업데이트 한다. λkφ 의 사후 분포를 예로 들면, π(λkφ ) π(λkφ ) Nk Y p(mi λkφ, λkτ, λk, µ1,, µnk ) i=1 이다. (supplemental material [Fu(2013)]) 2.3.3 위치 확률 행렬 P 의 추정 사후 위치 확률(P (N K))에 대한 예측 과정을 소개 하고자 한다. P 의 원소 pik = P r(zi = k M )는 i번째 객체가 군집 k에 분류될 확률을 뜻한다. Q 를 h번째 MCMC 표본 추출 과정에서의 사후 위치 확률 행렬이라고 한다. 그리고 K, α, α K α이면, γk G0, k = 1,, K (ω1,, ωk ) DirichletK (α,, α ), (2.25) P r(γi = γk ) = P r(zi = k) = ωk. h번째 MCMC 표본 추출과정에서의 디리크레 과정 혼합 모형은, (ω1,, ωk ) DirichletK h (α,, α ), P r(zi Mi {Zi = k w ) = ωk, (2.26) = k, Θk, Σk } NJR (Θk, Σk ). h번째 과정에서의 모수가 ξ 라고 하면, 사후 확률 행렬의 원소 qik 는 다음과 같은 분포에서 추출 할 수 있다. qik = P r(zi = k M, ξ ) (2.27) NJR (Mi Θk, Σk )ωk NJR (Mi Θk, Σk )DirichletK (ωk α,, α ). 15

지금까지 디리크레 과정 무작위 효과 혼합 모형의 이론 배경에 대해서 소개를 하였고, 3장에서는 시뮬레이션 데이터에 대한 분석을 살펴보도록 한다. 16

Chapter 3 시뮬레이션 자료분석 이 부분에서는 시뮬레이션 자료를 생성하고 MH 알고리즘을 이용하여 분석한 결과에 대해서 알아 보도록 한다. Table 3.1: K p λφ λτ 6 0.05 0.1 0.05 p λφ, λτ, λ, 각 군집의 평균과 분산 각 군집의 자료의 평균 각 군집의 자료의 분산 자료크기 0.5 0.1 0.4 20 0.1 0.5-1.0 0.1 20 0.05 0.5 0.2 2.5 0.4 20 0.05 0.5 0.2-1.5 0.4 20 0.05 0.3 0.3 0.5 0.2 20 0.05 0.3 0.3 0.1 1.0 20 λ 시뮬레이션 자료는 최대한 유전자 자료와 유사하게 생성하였다. 120개의 객 체를 18개의 시점(0,5,10,15,20,25,30,35,40,50,60,70,80,90,100,110,120,150)에서 17

4번 반복한 자료이며, 군집의 개수는 6개로 설정을 하였다. 즉, N = 120, J = 18, R = 4, K = 6인 자료이다. 각 군집의 자료들은 다중 정규 분포에서 추출하였고, 평균 벡터는 OU 과정을 통하여 추출하였다. 또한, OU 과정에서 필요하는 λφ, λτ, λ 의 참값은 각각 표3.1에서 제시 된 값을 사용하였다. 생성한 자료의 크기는 120 72이고, 6 객체 평균 벡터 행렬의 크기는 6 18이다. 생성된 시뮬레이션 자료는 그림3.1 과 같다. 여기서 각 그래프는 하나의 군집을 나타내고, 각 군집의 평균 벡터를 그린 그래프이다. Figure 3.1: 시뮬레이션 자료 시뮬레이션 자료를 디리크레 과정 무작위 효과 혼합 모형에 적용한 결과 총 8개의 군집으로 나눠졌고, 그 중 하나의 객체만 포함한 군집은 2개이다. 군집 분 석한 결과를 각 객체의 평균 벡터를 시각화하여 그림으로 나타났다(그림 3.2). 또한 각 자료의 군집의 참값과 비교하여 랜드 인덱스(Rand Index[Hubert and Arabic (1985)])를 계산한 결과 0.914(1에 가까워 질 수록 좋다)로 비교적 정확한 결과를 얻었다고 볼 수 있다. 단, 시간이 다소 오래 걸린 다는 것을 단점으로 볼 수 있으며, 이번 시뮬레이션 자료의 군집 분석 과정에서 10,000번의 반복 시행에 18

소요되는 시간은 약 6시간이었다. Figure 3.2: 시뮬레이션 자료 군집 분석 결과(평균 벡터) 그림 3.2에서 검은 색이 아닌 라인 들은 각 군집 평균 벡터의 사후 분포 평 균의 추정치이다. 또한, 현재 각 군집에 분류된 객체들의 평균 백터는 분포의 일부분이므로, 이들의 사후 분포의 평균이 꼭 객체를 나타나는 라인들의 중심에 위치하지 않아도 된다. 그림 3.3에서는 각 군집의 평균 벡터를 추출하는 과정에 p 서 사용되는 λφ, λτ, λ 의 추정 값이다. 19

Figure 3.3: 시뮬레이션자료군집분석결과 ( λ φ, λ τ, λ ɛ ) 20

Chapter 4 유전자 자료분석 여기서 다루게 되는 유전자 자료는 Housden(2011)에서 소개한 자료(전체 자료 중의 일부 만 사용)이다. 총 18 4 = 72개의 초파리 성체 근육세포 마이 크로 어레이 유전자 자료를 로그 변환 시킨 값이다. 또한 자료의 반복 측정은 4개의 생물학 복제물을 통하여 이루어졌으며, 163개의 차별적 유전자에 대해 군집 분석을 하였다. Figure 4.1: 유전자 자료 군집 분석 결과(평균 벡터) 유전자 자료를 군집 분석한 결과, 총 15개의 군집으로 나누어졌고, 그 중 21

3개의 군집은 하나의 유전자 만 포함되어 있다. 그림 4.1에서는 2개 이상의 유 전자를 포함하는 군집에 대한 그림이다. 군집 12에서는 군집 내 유전자의 평균 벡터의 사후 분포(다변량 정규 분포)에서 평균(Θk )의 추정값의 그래프는 평균 벡터 들의 중심에 위치하지 않았다. 이는 위치 확률에 의하여, 이 분포를 따르는 여기에서 나타나지 않는 자료들을 많이 사용했기 때문이다. p Figure 4.2: 유전자 자료 군집 분석 결과( λφ, λτ, λ ) 그림 4.2에서는 3개 표준 편차 p λφ, λτ, λ 에 대한 추정 결과를 나타나는 그래프이다. 3개의 표준편차 모두 각 군집 사이에서 큰 차이를 띄고 있으며, 대 부분의 군집 안의 무작위 효과를 나타나는 표준 편차의 값이 작다는 것을 알 수 있다. 그림 4.3 은 사후 위치 확률 행렬에 대한 PCA를 기반으로 군집 분석한 결과를 시각화한 결과이다. 또한 그림 4.3에서 각 군집을 나타나는 색은 그림 4.1, 그림 4.2와 같다. 확률 행렬의 각 행은 유전자 객체의 군집 위치에 대한 분포이기 때 문에 PCA그림(그림 4.3)은 각 유전자 객체의 현재 위치와 다른 군집의 관계를 나타나는 그림이기도 한다. 비슷한 사후 확률을 가진 유전자들은 가까이 위치해 22

Figure 4.3: 유전자자료군집분석결과 (PCA) 있고, 확률값이 0.8 이상인유전자들은하나의군집에분류되어있다. 하지만 2 개의군집에동시분류된유전자들은군집사이의경계에위치하고있음을알 수있다. 23

Chapter 5 결론 본 논문에서는 우선 모형에 사용되는 디리크레 과정에 대한 이론 배경을 간략 하게 소개하였고, 반복 측정한 마이크로 어레이 유전자 자료를 분석하는 디리크 레 과정 무작의 효과 혼합 모형(Fu (2013))을 소개하였다. 그리고 Housden(2011) 에서 소개한 초파리 성체 근육세포 유전자 자료의 일부에 적용하여 분석하였다. 군집의 개수가 명확하지 않는 유전자 자료의 분석에서는 베이지안 비모수 방 법론이 널리 사용되고 있으며, 유전자 분석을 통한 질병의 진단 및 예측에 큰 기대를 할 수 있다. 모형에서는 시간의 간격에 대한 제한이 없으므로 적용 할 수 있는 데이터 유형이 국한 되어 있지 않는 점에서, 군집의 개수가 명확하지 않고, 여러 시점에서 반복으로 측정된 자료를 분석할 수 있는 장점이 있다. 단점 은 MCMC알고리즘을 사용함으로서 분석 시간이 오래 걸린다는 점이다. 단점을 보완하여 사용하면 정확하고 빠른 군집 분석 및 예측을 할 수 있을 것이다. 24

참고문헌 J.K. Ghosh, and R.V.Ramamoorthi(2003), Bayesian Nonparametrics, Springer, 87 116. Qiuyan Fu(2013), Baysian Clustering of Replicated time-course gene expression data with Weak Signals, Cambridge. John Paisley, (2010), A Simple Proof of the Stick-Breaking Construction of the Dirichlet Process, Princeton. Yee Whye Teh, Dirichlet Process, University College London. Emin Orhan (2012), Dirichlet Processes, Rochester education. Neal (2000), Markov Chain Sampling Methods for Dirichlet Process Mixture Models, Journal of Computational and Graphical Statistics, 249 265. 25