The Korean Journal of Applied Statistics 2021, Vol. 34, No. 1, 9 23 DOI: http://dx.doi.org/10.5351/kjas.2021.34.1.009 A variational Bayes method for pharmacokinetic model Sun Parka, Seongil Jo 1,b, Woojoo Leec a Department of Statistics, Jeonbuk National University; b Department of Statistics, Inha University; c Graduate School of Public Health, Seoul National University Abstract In the following paper we introduce a variational Bayes method that approximates posterior distributions with mean-field method. In particular, we introduce automatic differentiation variation inference (ADVI), which approximates joint posterior distributions using the product of Gaussian distributions after transforming parameters into real coordinate space, and then apply it to pharmacokinetic models that are models for the study of the time course of drug absorption, distribution, metabolism and excretion. We analyze real data sets using ADVI and compare the results with those based on Markov chain Monte Carlo. We implement the algorithms using Stan. Keywords: automatic differentiation variational inference, markov chain monte carlo, pharmacokinetic models, stan, variational bayes 1. 서론 베이지안 추론(Bayesian inference)은 자료를 분석하기 위한 방법의 일환으로 관측치가 주어졌을 때 관심있는 모수의 조건부 분포인 사후 분포(posterior distribution)를 계산하여 추론하는 방법이다 (Murphy, 2012; Gelman 등, 2013). 베이지안 추론의 장점 중 하나는 모수의 추론에서 확률분포를 이용하기 때문에 자연스럽게 불확 실성(uncertainty)를 고려할 수 있다는 것이다. 그러나 모형이 복잡하거나 모수의 차원이 높아지게 되면 사후 분포의 계산이 어려워져 직접 구할 수 없거나, 구하더라도 잘 알려진 형태로 얻어지지 않는다는 것이 알려져 있다 (Bishop, 2006). 따라서 실제 연구에서는 사후 분포에 비례하는 확률분포에서 표본을 생성시켜 추론하 는 방법인 마코프 체인 몬테 카를로(Markov chain Monte Carlo; MCMC) 방법을 사용해왔다 (Hastings, 1970; Gelfand와 Smith, 1990; Robert와 Casella, 2013). MCMC는 비교적 정확한 결과를 내지만 막대한 양의 계산식이 수반되기 때문에 자료의 양이 많아지거나 모형이 복잡해지면 계산 시간이 매우 길어지게 된다는 단점을 가지고 있다 (Jordan 등, 1999). 따라서 빠른 결과와 더불어 신속한 의사 결정이 필요한 연구에서는 MCMC를 이용한 방법이 적절하지 않을 수 있다. 최근 MCMC가 가진 이러한 문제점을 해결하기 위해 변분 방법(variational method)을 활용한 베이지안 추론인 변분 추론(variational inference)이 활발히 개발 되고 있다. 변분 추론은 Jordan 등 (1999)에 의해 제안된 방법으로 계산이 쉬운 분포 중 사후 분포에 가장 근접하는 특정 분포를 찾아 추론하는 방법이다. 이 방법의 가장 큰 장점은 계산이 어려운 사후 분포를 대신해 계산이 This paper is based on part of Sun Park s master thesis. Research of Seongil Jo was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (NRF-2020R1C1C1A01013338). 1 Corresponding author: Department of Statistics, Inha University, 100 Inha-ro Michuhol-gu, Incheon 22212, Korea. E-mail: statjs@inha.ac.kr
Sun Park, Seongil Jo, Woojoo Lee 10 비교적 용이한 근사 분포를 이용하기 때문에 계산 소요시간이 단축되어 고차원의 복잡한 모형 뿐 아니라 대용 량의 자료에서도 빠르게 사용할 수 있다는 점이다 (Wainwright과 Jordan, 2008; Blei와 Jordan, 2004). 하지만 사후 분포에 근접하는 분포를 찾는 일은 쉽지 않다. 최근 이를 해결하기 위해 많은 연구가 진행되고 있으며 Kucukelbir 등 (2015)은 모수를 실수 공간으로 변환하여 사후 분포를 가우시안 분포로 근사하는 방법인 자 동 미분 변분 추론(automatic differentiation variational inference; ADVI)을 제안하였다. 따라서 본 논문에서는 ADVI에 대해 살펴보고 약물동태학 모형(pharmacokinetic models)에 적용하여 실제 문제에 있어서 MCMC의 대안으로 사용이 가능한지 살펴보고자 한다. 약물동태학은 환자에게 약물을 투여한 후 체내에서의 약물의 움직임을 연구로 환자의 약물 투여 계획을 세울 때 이용된다. 약물은 환자 개개인의 특성에 따라 다르게 작용되기 때문에 약물이 안전하면서 가장 효율 적으로 작용할 수 있도록 투여량, 투여 간격 등에 대한 계획을 세워야 한다. 여러 요인들이 체내에서 약물의 움직임에 영향을 주어 약물의 효과에 변화를 가져올 수 있기 때문에 상황에 따라 환자의 상태를 살피면서 다음 투여가 언제 이루어져야 하는지 빠르게 예측해야하는 경우가 발생한다. 이런 경우, 약동학에서는 사전 정보 를 이용하여 비교적 빠른 분석이 가능한 베이지안 추론을 많이 사용한다. 그러나 MCMC를 사용하면 계산은 정확할지라도 모형의 복잡성, 자료의 양에 따라 계산에 소요되는 시간이 길어지게 되기 때문에 투여 계획을 세우는데 문제가 생길 수 있다. 따라서 본 연구를 통해 ADVI의 성능이 입증된다면 실제 상황에서 빠르게 약물 투여 계획을 세울 수 있어 환자 치료에 상당한 도움이 될 것이다. 본 논문의 나머지 부분은 다음과 같이 구성된다. 먼저 2절에서는 변분 방법을 활용한 베이지안 추론인 변분 추론의 기본적인 개념과 변분 베이즈 방법, 그리고 자동 미분 변분 추론에 대해 설명하고, 3절에서 약 물동태학 모형에 적용을 통해 MCMC방법과의 비교 결과를 제시한다. 그 다음 마지막으로 4절에서 ADVI의 장점과 단점에 대해 기술하고 논문을 마친다. 2. 방법론 본 절에서는 변분 추론과 변분 베이즈 방법, 그리고 자동 미분 변분 추론에 대해 간략히 설명하고 이를 구현할 수 있는 Stan (Stan Development Team, 2017)을 소개한다. 2.1. 변분 추론 변분 추론이란 계산이 어려운 사후 분포를 다루기 쉬운 분포로 근사하여 추론하는 방법으로 최근 활발히 연 구가 되고 있다 (Jordan 등, 1999; Blei와 Jordan, 2004). 이를 구체적으로 설명하기 위해 먼저 관측된 변수를 x, 관심이 있는 모수를 θ Θ라 하고 사후 분포를 p(θ x)라 하자. 그때, 변분 추론은 아래와 같이 쿨백-라이블러 발산(Kullback-Leibler divergence; KL divergence)을 최소화 하는 분포 q (θ)를 찾아 사후 분포의 근사 분포로 하여 추론을 하는 것이다. q (θ) = arg min KL (q(θ) k p(θ x)), (2.1) q(θ) Q 여기서 Q = {q(θ); θ Θ}는 근사 분포의 집합을 나타내고, KL ( q(θ) k p(θ x) ) = Z q(θ) log q(θ) dθ. p(θ x) (2.2) 식 (2.1)을 최소화하는 q (θ)를 찾는 것은 개념적으로 쉽다. 그러나 위의 쿨백-라이블러 발산식 (2.2)을
A variational Bayes method for pharmacokinetic model 11 정리해보면, Z q(θ) dθ p(θ x) Z Z p(θ, x) = q(θ) log q(θ) dθ q(θ) log dθ p(x) Z Z = q(θ) log q(θ) dθ q(θ) log p(θ, x) dθ + log p(x) KL ( q(θ) k p(θ x) ) = q(θ) log = E log q(θ) E log p(θ, x) + log p(x) 으로 로그를 취한 정규화 상수(normalizing constant) log p(x)가 포함되어있기 때문에 계산이 어렵다. 따라서 실제 적용에 있어서는 아래와 같이 쿨백-라이블러 발산식에 다시 한번 정규화 상수를 더한 Evidence Lower Bound (ELBO)를 최대화 하는 근사 분포를 찾는다. L (q(θ)) = E log q(θ) E log p(θ, x) + log p(x) + log p(x) = E log p(θ, x) E log q(θ). (2.3) 위 식에서 기대값은 q(θ)에 대한 것이다. 변분 추론의 좀 더 자세한 내용에 대해서는 Wainwright과 Jordan (2008)를 참고하기 바란다. 2.2. 변분 베이즈 변분 베이즈(variational Bayes; VB)방법은 평균장 근사(mean-field approximation)를 이용하여 사후 분포를 근 사하는 방법 (Ghahramani와 Beal, 2001; Beal, 2003)으로 관측치 x에 대한 모수가 p차원의 벡터 θ = (θ1,..., θ p )> 라 할 때, 근사 분포 집합 Q의 원소를 각 모수 θ j, j = 1,..., p에 대한 확률 분포 q j (θ j )의 곱으로 가정한다. 즉, 다음과 같이 인수분해 형태(factorized form)로 정의하여 사용한다. q(θ) = p Y q j (θ j ), j = 1,..., p. (2.4) j=1 참고로 위 식 (2.4)는 모수에 대한 확률 분포들이 서로 독립임을 의미한다. 변분 베이즈 방법은 독립 이외의 추가적인 가정은 필요없으며 q j (θ j )에 대해서도 별다른 제약이 필요없다 (Bishop, 2006). 또한 모수에 대한 확률 분포들이 서로 독립임을 가정하기 때문에 계산식이 쉬워져 속도가 빠르고 대용량의 자료 분석에 장점을 갖는다. 다만, 모수들 간의 상관성을 고려하지 않기 때문에 불확실성 (uncertainty)이 과소 추정(under-estimation)된다는 문제점이 있다 (Bishop, 2016). 2.3. 자동 미분 변분 추론 변분 추론을 위해서는 모수들에 대해 근사 분포를 특정해야 하는데 그 과정이 쉽지 않다. 이를 해결하기 위해 Kucukelbir 등 (2015, 2017)는 ADVI방법을 제안하였다. ADVI의 기본적인 아이디어는 실수 공간에서 정의 되지 않은 모수들을 실수 공간의 지지집합(support)을 가지도록 일대일 변환 후에 변환된 모수의 사후 분포 를 가우시안 분포(Gaussian distribution)로 근사하는 것이다. 이때, ELBO의 최대화는 자동 미분(automatic differentiation)과 확률적 최적화(stochastic optimization)를 결합하여 실시한다. ADVI 방법 중 평균장 방법을 기반으로 하는 평균장 가우시안(mean-field Gaussian)방법에 대해 구체적으로 나타내면 다음과 같다.
Sun Park, Seongil Jo, Woojoo Lee 12 1. 먼저, 모수의 지지집합(support)이 실수 공간에 존재하도록 아래와 같이 일대일 미분 가능 함수를 정의한 후 T : supp (p(θ)) R p, (2.5) 변수 변환을 실시한다. 즉, 모수가 p차원의 벡터 θ = (θ1,..., θ p )T 일 때, 식 (2.5)에서 정의된 T ( )를 이용하여 ζ = T (θ)로 변환한다. 이때, 변환된 모수의 결합 밀도 함수는 아래와 같이 표현된다. p(x, ζ) = p x, T 1 (ζ) det JT 1 (ζ). (2.6) 위 식에서 JT 1 (ζ)는 T 1 의 자코비안(Jacobian)을 나타낸다. 참고로 자코비안은 변환된 확률 밀도 함수가 적분해서 1이 되는 성질을 만족시킬 수 있게 한다 (Olive, 2014). 2. 두 번째 단계에서는 변환 된 모수 ζ의 사후 분포를 가우시안 분포로 근사한다. 구체적으로, p (ζ x)의 근사 분포로 아래와 같이 인수분해 형태의 분포를 고려한 후, q(ζ; φ) = p Y N ζ j µ j, σ2j, (2.7) j=1 여기서 φ = (µ1,..., µ p, σ21,..., σ2p )T 이다. 다음과 같이 주어지는 ELBO를 최대화 하는 평균과 분산을 찾는 다. p h i p X 1 + log(2π) + log σ j. L(φ) = E log p x, T 1 (ζ) + log det JT 1 (ζ) + 2 j=1 (2.8) 3. ELBO 식 (2.8)를 최대화하는 계산식은 적분이 복잡하고 분산 σ2 이 항상 0보다 크다는 제약조건 때문에 계산이 어렵다. 따라서 Kucukelbir 등 (2017)는 추가적으로 elliptical standardization (Ha rdle과 Simar, 2012) 을 실시하여 그 문제를 해결하였다. 즉, 제약조건이 있는 분산에 대해 로그 변환 ω j = log(σ j ), j = 1,..., p 후 표준화 η = S µ,ω (ζ) = diag(exp(ω) 1 )(ζ µ)를 통해 평균이 0 분산이 1로 표준화된 가우시안 분포를 활 용하여 최적화를 실시하였다. 이러한 방법은 확률적 최적화를 위한 기울기(gradient) 계산을 쉽게 만들기 때문에 변분 추론에서 매우 유용한 방법이다. Elliptical standardization을 통한 구체적인 근사 분포는 q(η; 0, I) = N(η 0, I) = p Y N(η j 0, 1) j=1 이고 ELBO는 다음과 같이 주어진다. p X 1 1 L (µ, ω) = EN(η 0,I) log p x, T 1 S µ,ω (η) + log det JT 1 S µ,ω (η) + ω j. j=1 지금까지 설명한 평균장 가우시안 방법 이외에 모수들의 상관성을 고려할 수 있는 완전 계수 가우시안 (full-rank Gaussian) 분포를 이용하는 방법과 ADVI의 이론적인 성질에 대해서는 Kucukelbir 등 (2017)를 참고 하기 바란다. 2.4. Stan 본 논문에서는 소개된 ADVI를 구현하기 위해 Stan (Stan Development Team, 2017)을 사용하였다. Stan은 베이지안 분석에서 통계 모형을 구현하기 위한 C++ 기반의 확률적 프로그래밍 언어로 통계 모형 자체를 코 드로 작성하기 때문에 복잡한 통계 모형도 쉽게 구현이 가능하다는 장점이 있다. 또한 Stan은 MCMC 방법 중
A variational Bayes method for pharmacokinetic model 13 Figure 1: One-compartment model for oral administration. 해밀토니안 몬테 칼로(Hamiltonian Monte Calro) 기반의 No-U-Turn Sampler (Hoffman과 Gelman, 2014)를 사 용한 베이지안 추론도 가능하여 최근 널리 사용되고 있는 프로그램 중 하나이다. Stan은 R에서 관련 패키지인 RStan을 설치 한 후 사용할 수 있다. 참고로 R에서 Stan을 이용 할 경우 결과를 표나 그래프를 통해 시각적으로 쉽게 이해할 수 있으며 새로 추출한 표본을 자유롭게 사용할 수 있다는 장점이 있다 (Matcuura, 2016). 3. 분석 3.1. 약물동태학 약물동태학이란 체내에 투여된 약물이 시간의 경과에 따라 체내에서 흡수(absorption), 분포(distribution), 대 사(metabolism) 및 제거(excretion)되기까지 일련의 과정을 동태학적 입장에서 이해하고자하는 학문으로 임상 시험 및 약물을 통한 환자 치료에 중요하게 사용되고 있다 (Shim, 1994; Wakefield, 2013). 체내에서 약물의 움직임을 예측하면 신약 개발 시 임상 시험을 통해 안전하면서도 효율적인 신약을 개발 할 수 있을 뿐만 아니라 환자에게 투여된 약물이 안전하면서 효율적으로 작용하도록 적절한 투여량과 투여 간격 및 투여 제형에 대한 계획을 세울 때에도 도움을 준다 (Shim, 1994; Choi, 2000). 특히 약동학 실험을 통해 투여 계획을 세울 때에는 성별, 키, 몸무게 등을 포함한 환자 개개인의 특성을 살피어 지속적인 모니터 링을 통해 빠르게 예측해야하기 때문에 사전 정보를 이용한 베이지안 추론이 많이 사용되고 있다 (Wakefield, 2013). 약물동태학에서는 약물의 체내 이행 과정을 설명하기 위하여 농도 평형 관계가 존재하는 조직들을 하나로 묶어 구획(compartment)이라 한다 (Shim, 1994). 아래에 주어진 Figure 1은 투약 후 체내에서 약물의 움직임을 가장 간단하게 설명하는 일구획 모형(1-compartment model)을 도식화한 그림이다 (Godfrey, 1983). 일구획 모형은 약물이 주입되는 구획 0과 약물이 흡수되는 구획 1로 이루어져있다. 여기서 구획 1은 혈장으로 blood compartment 라고도 하며 혈장으로 흡수된 약물이 분해가 되어 체내에서 소실되는 과정을 거친다. 따라서 채취된 혈액 표본의 혈장에서는 구획 1에서 발생하는 약물의 분해 과정을 확인할 수 있다 (Wakefield, 2013). 일구획 모형에서 ωk (t)는 t시점에서의 구획 k (= 0, 1)내 약물의 총량이라 할 때, 두 구획간의 약물의 흐름은 총량 ωk (t)를 t로 미분함으로써 구할 수 있다. 구체적으로, 단위시간 단위부피 당 약물이 구획 0에서 구획 1로 ka > 0만큼 흘러들어가고 단위시간 단위부피 당 약물이 구획 1(혈장)내에서 ke > 0만큼 분해된다고 할 때, 미분방정식은 d ω0 = ka ω0, dt d ω1 = ka ω0 ke ω1 dt (3.1) (3.2)
Sun Park, Seongil Jo, Woojoo Lee 14 Table 1: Concentraion of the drug theophylline as a function of time, obtained from a patient who was administered an oral dose of size 4.53 mg/kg Observation number i 1 2 3 4 5 6 7 8 9 10 Time (hours) xi 0.27 0.58 1.02 2.02 3.62 5.08 7.07 9.00 12.15 24.17 Concentration (mg/liter) yi 4.4 6.9 8.2 7.8 7.5 6.2 5.3 4.9 3.7 1.05 로 주어지고, t = 0일 때 초기 복용량 ω0 (0)을 D로 정하고 미분방정식 (3.1)와 (3.2)를 풀어 x시점에서의 혈중 약물 총량 ω1 (x)을 구한 후 구획 1의 부피인 V로 나눔으로써 표준화하여 혈중 약물 농도 µ(x)를 구하면 아래와 같이 주어진다: µ(x) = Dka exp( ke x) exp( ka x). V(ka ke ) 위 식으로부터 혈중 약물 농도를 측정함으로써 혈액에서 약물이 분해되는 속도인 소실속도 상수(elimination rate constant) ke 를 구할 수 있다. 그러나 소실속도 상수는 어디까지나 속도 상수일뿐이지 그 자체로 속도를 나타내는 것이 아니기 때문에 그 의미를 해석하기에 어려움이 있다. 따라서 혈중 약물 농도가 절반 이 되는 시간을 나타내는 소실 반감기(half-life)와 부피를 고려한 소실의 정도를 나타내는 청소율(clearence) 를 사용하여 약물이 체내로부터 소실되는 속도의 척도로 삼는다. 소실 반감기와 청소율은 각각 t1/2 과 CL로 나타내며 구하는 식은 아래와 같다. t1 = 2 log 2, ke CL = V ke. 3.2. 1명의 환자에 대한 경시적 자료 본 절에서는 실제 자료를 가지고 MCMC와 ADVI 두 방법으로 약물동태학 모형에서의 모수를 추정하고 비교 한다. 자료는 체중이 70.5kg인 환자 1명에 대한 항정신병약물인 테오필린(Theophylline)을 복용한 후 10개의 시점에서 혈중 약물 농도를 측정한 자료로 Wakefield (2013)에서 가져왔으며 Upton 등 (1982)에 의해 연구 된 자료이다. 테오필린의 초기 복용량은 kg 당 4.53mg으로 대상자의 경우 70.5kg 4.53mg = 319.365mg을 복용하였다. 시간의 흐름에 따른 약물의 농도에 대한 자료의 구체적인 값은 Table 1에 나타나있다. 3.2.1. 모형 xi 를 혈중 약물 농도를 측정한 시점이라 하고 yi 를 xi 시점에서 측정한 혈중 약물의 농도라 하자. 그때, 환자가 약 물을 복용(D) 한 후 시간의 경과에 따라 약물의 혈중 농도는 아래와 같이 비선형 회귀모형(nonlinear regression
A variational Bayes method for pharmacokinetic model 15 Table 2: Comparison results for analyzing one-compartment model using Stan Parameters σ ka ke CL V t1/2 소요시간(sec) MCMC 0.26 2.46 2.78 34.30 8.55 ADVI 0.20 0.00 0.13 0.93 0.55 0.26 2.48 2.77 33.94 8.48 183.37 0.13 0.00 0.04 0.81 0.15 2.38 model)을 고려하고 yi N µi, σ2, µi = D ka V (ka ke ) i = 1,..., 10, e ke xi e ka xi, i = 1,..., 10, 각 모수에 대한 사전 분포는 아래와 같이 설정하였다. σ2 Uni f (0.01, 100), log(ka ) N( 2.52, 100), log(ke ) N(0.4, 100), log(cl) N( 3.25, 100). 위의 사전분포에서 ka 와 ke 는 속도 상수이기 때문에 로그를 취한 log(ka )과 log(ke )이 가우시안 분포를 따른 다고 설정하여 항상 양수가 되도록 하였으며, CL 역시 약물이 체내에서 소실되는 속도를 나타내기 때문에 항상 양수이므로 로그를 취한 log(cl)이 가우시안 분포를 따른다고 설정하였다. 초모수(hyperparameter)들의 값은 Wakefield (2013)에서 일반화 추정 방정식(generalized estimating equations; GEE)으로 계산한 모수의 추정치를 참고하였다. 3.2.2. 분석 일구획 모형의 모수 ka, ke, CL, V, t1/2 를 추정하기 위하여 Stan을 이용하였다. 이때, MCMC의 경우 서로 다 른 초기값을 가진 4개의 chain을 사용하였으며, 각 chain마다 총 14,000개의 표본을 추출하여 초기 7,000개의 표본을 버리고 남은 7,000개의 표본을 분석에 사용하였다. 수렴여부에 대해서는 Gelman과 Rubin (1992)에 의 해 제안된 potential scale reduction statistic값을 이용하였다. ADVI의 경우에는 140,000번의 반복을 하였으며 ELBO 값의 변화를 이용하여 수렴여부를 판단하였다. 참고로 MCMC의 경우 potential scale reduction statistic 값이 1에 가까우면 수렴했다고 가정하고 ADVI의 경우에는 ELBO 값의 변화가 거의 없으면 수렴했다고 가정 한다. Table 2는 MCMC와 ADVI 두 가지 방법으로 모수를 추정한 결과 비교표이다. 추정된 모수들의 기댓값을 확인해보면, 두 방법으로 부터 계산된 사후 분포의 평균이 비슷하게 추정되었음을 알 수 있다. 그러나 수렴까 지의 계산 소요 시간이 MCMC는 183.37초, ADVI는 2.38초로 약 76.7배 가량 계산 소요 시간이 단축되었음을 확인할 수 있다. Figure 2는 각 모수들의 사후 확률 분포를 나타낸다. 파란선이 MCMC로 계산된 사후 분포, 빨간선이 ADVI 로 근사한 사후 분포이다. 그림에서 볼 수 있듯이 사후 분포의 평균은 비슷하게 추정되고 있다. 그러나 MCMC
Sun Park, Seongil Jo, Woojoo Lee 16 (a) σ (b) ka (c) ke (d) CL (e) V (f) t1/2 Figure 2: Posterior distributions for parameters: blue lines represent posteriors from MCMC, red lines denote posteriors from ADVI. Figure 3: Concentration and 95% credible interval estimated as a function of time. Dots represent real observa- tions. 보다 ADVI의 분산이 더 적게 추정됨을 알 수 있다. Figure 3은 실제 혈중 약물 농도(concentration) y의 평균 반응에 대한 추정치와 95% 신용구간(credible interval)을 나타낸 것이다. Figure 2에서와 마찬가지로 파란선이 MCMC로 추정한 결과이고, 빨간선이 ADVI 로 추정한 결과를 나타낸다. 점은 실제 혈중 약물 농도 y를 나타낸다. 그림에서 볼 수 있듯이 두 방법 모두 y 가 평균 반응을 잘 추정했다고 할 수 있다. 그러나 95% 신용구간을 살펴보면 MCMC보다 ADVI가 신뢰 구 간이 더 좁게 추정됨을 확인할 수가 있다. 이는 ADVI가 MCMC에 비해 불확실성(uncertainty)을 과소 추정
A variational Bayes method for pharmacokinetic model 17 Table 3: The oral dose(mg) of theophylline for each subject Subject ID Initial dose 1 320.00 2 320.10 3 319.80 4 320.65 5 318.56 6 319.37 Subject ID Initial dose 7 319.88 8 319.96 9 320.00 10 319.77 11 319.36 12 267.84 (under-estimation)하는 경향이 있다는 것을 의미한다. 3.3. 12명의 환자에 대한 경시적 자료 본 절에서는 12명의 환자에 대해서 테오필린을 복용한 후 시간의 경과에 따른 혈중 약물 농도를 측정한 자료를 이용하여 비교분석을 실시한다. 초기 복용량은 Table 3에 나타나있다. 3.3.1. 모형 i번째 환자가 복용한 약물의 양을 Di, j번째 시점 xi j 에서 측정한 혈중 약물의 농도를 yi j 라 하자. 그때 다음의 계층적 모형(hierarchical model)을 통해 약물의 혈중 농도를 모형화 하고 yi j N µi j, σ2, i = 1,..., 12, j = 1,... 10, Di kai µi j = e kei xi j e kai xi j, i = 1,..., 12, j = 1,... 10, Vi (kai kei ) CLi, i = 1,..., 12, Vi = kei log 2 t1 =, i = 1,..., 12, 2i kei 사전 분포를 다음과 같이 계층적으로 설정하였다. log(kai ) N lka, σ2ka, i = 1,..., 12, log(kei ) N lke, σ2ke, i = 1,..., 12, log(cli ) N lcl, σ2cl, i = 1,..., 12, lka Uni f ( 3.5, 0.5), lke Uni f ( 2, 2.5), lcl Uni f ( 3, 3), σ2 Uni f (0, 100), σ2ke Uni f (0, 100), σ2ka Uni f (0, 100), σ2cl Uni f (0, 100). 3.3.2. 분석 Table 4와 Table 5는 MCMC와 ADVI 두 가지 방법으로 모형을 구현한 결과 비교표이다. Table 4는 12명 환자 전체에 대한 모수의 추정치를 나타내고, Table 5는 환자 개개인에 대한 모수의 추정치를 나타낸다. 두 표에서
Sun Park, Seongil Jo, Woojoo Lee 18 Table 4: Comparison results obtained from MCMC and ADVI Parameters σ lka lke lcl σka σke σcl 소요시간(sec) MCMC 0.71 0.48 2.46 1.00 0.69 0.03 0.22 ADVI 0.05 0.21 0.05 0.07 0.17 0.03 0.05 0.75 0.54 2.46 1.00 0.63 0.04 0.22 994.77 0.05 0.18 0.01 0.05 0.14 0.01 0.05 7.57 Table 5: Comparison results obtained from MCMC and ADVI for each patient Parameters CL1 CL2 CL3 CL4 CL5 CL6 CL7 CL8 CL9 CL10 CL11 CL12 t1/21 t1/22 t1/23 t1/24 t1/25 t1/26 t1/27 t1/28 t1/29 t1/210 t1/211 t1/212 MCMC 2.17 2.91 2.86 2.69 2.34 3.63 3.04 3.11 3.26 2.07 3.32 1.94 8.33 8.14 8.21 8.19 8.16 8.22 8.20 8.20 8.21 8.21 8.18 8.09 ADVI 0.12 0.18 0.18 0.16 0.13 0.26 0.20 0.20 0.21 0.11 0.22 0.11 0.58 0.54 0.53 0.54 0.52 0.55 0.54 0.54 0.53 0.53 0.54 0.55 2.16 2.90 2.85 2.68 2.32 3.59 3.04 3.09 3.26 2.05 3.28 1.93 8.35 8.18 8.22 8.18 8.19 8.24 8.23 8.20 8.28 8.15 8.18 8.13 0.13 0.12 0.11 0.21 0.16 0.15 0.14 0.07 0.15 0.07 0.24 0.25 0.26 0.27 0.24 0.29 0.29 0.28 0.27 0.25 0.27 0.26 볼 수 있듯이 추정된 모수들의 기댓값을 확인해보면, 두 방법이 비슷하게 추정되었음을 알 수 있다. 그러나 계산 소요 시간이 MCMC는 994.77초, ADVI는 7.57초로 약 131.4배 가량 계산 소요 시간이 단축되었음을 확인할 수 있다. Figure 4와 Figure 5는 각 모수들의 사후 확률 분포를 보여준다. 파란선이 MCMC으로 추정한 사후 분 포, 빨간선이 ADVI로 추정한 사후 분포이다. 그림에서 볼수 있듯이 사후 분포의 평균은 비슷하게 추정되나 MCMC보다 ADVI의 분산이 더 작게 추정 됨을 알 수 있다.
A variational Bayes method for pharmacokinetic model 19 Figure 4: Posterior distributions of CL for 12 patients. Figure 5: Posterior distributions of t1/2 for 12 patients. Figure 6은 주어진 실제 혈중 약물 농도와 MCMC와 ADVI로 추정한 평균 반응 함수, 95% 신용구간를 나타낸다. 기존과 마찬가지로 파란선이 MCMC, 빨간선이 ADVI, 점은 실제 혈중 약물 농도 y를 나타낸다. 그림에서 보이는대로 두 방법 모두 혈중 약물 농도를 평균적으로 잘 추정했다고 할 수 있다. 그러나 95% 신용구간을 살펴보면 MCMC보다 ADVI가 신뢰 구간이 더 좁게 추정됨을 확인할 수가 있다. 이는 ADVI가 불확실성을 과소 추정 하는 경향이 있다는 것을 의미한다.
Sun Park, Seongil Jo, Woojoo Lee 20 Figure 6: Concentration and 95% credible interval estimated as a function of time for each patient. Dots represent real observations. Table 6: Comparison results for selecting prior distribution Parameters σ ka ke CL V t1/2 Informative prior 0.26 2.48 0.13 0.00 2.77 0.04 33.94 0.81 8.48 0.15 Non-informative prior 0.28 0.07 2.48 0.14 0.00 2.79 0.04 34.35 0.83 8.52 0.17 3.4. 민감도 분석 본 절에서는 추가적으로 사전분포가 모형의 추정 결과에 어떤 영향을 미치는지 확인하기 위한 민감도 분 석을 실시한다. 이를 위해 자료와 모형은 앞서 설명한 1명의 환자에 대한 것을 사용하고, 사전분포로는 본 논문에서 사용한 초모수를 갖는 정보적 사전분포(informative prior)와 초모수를 갖지 않는 무정보적 사전분포 (non-informative prior)를 고려한다. Table 6은 민감도 분석에 대한 결과를 나타낸다. 표로부터 모수 추정 결과를 비교해보면 그 값이 비슷하 다는 것을 알 수 있다. 이는 해당 모형의 추정 결과가 사전분포의 선택에 민감하게 반응하지 않는다는 것을 의미한다. 4. 결론 본 논문에서는 사후 분포의 계산 속도를 향상 시킬수 있는 자동 미분 변분 추론 ADVI에 대해 살펴보고 이 를 약물동태학 모형 중 1-컴파트먼트 모형에 적용하여 MCMC와의 성능을 비교하였다. 두 방법으로 구현한
A variational Bayes method for pharmacokinetic model 21 모형 모두 환자의 시간에 따라 측정한 혈중 약물의 농도를 평균적으로 잘 예측하였다고 판단 할 수 있었다. 그러나 결과를 도출하기까지의 계산 속도 즉, 사후 분포에 수렴하기 까지의 계산 속도에서 ADVI가 MCMC 에 비해 1명인 경우 약 76.7배, 12명인 경우 약 131.4배 가량 단축되었다. 이는 모형이 복잡해지거나 자료의 양이 많아질수록 MCMC와 ADVI의 계산 속도 차이가 더 커진다는 것을 나타낸다. 따라서 빠른 결과를 필요로 하는 분석 상황에서는 MCMC보다 ADVI가 더 적절하다고 판단된다. 다만, 평균장 가우시안 방법을 기반으로 하는 ADVI를 이용하여 사후 분포를 근사할 경우 모수들의 상관성을 고려하지 않기 때문에 불확실성이 과소 추정되는 문제점이 있음을 확인하였다. 따라서 불확실성을 중요하게 고려해야 하는 문제에서는 평균장 가우 시안 방법을 기반으로 하기 보다는 다소 시간이 걸리더라도 다른 방법 예를 들어 완전 계수 가우시안 ADVI (Kucukelbir 등, 2015, 2017)를 사용할 것을 추천한다. References Beal, M. J. (2003). Variational algorithms for approximate Bayesian inference, Ph. D. thesis, University College London, UK. Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Springer. Blei, D. M. and Jordan, M. I. (2004). Variational methods for the dirichlet processes. In Proceedings of the twenty-first international conference on Machine learning, pp. 12. Choi, J. S. (2000). Application of Pharmacokinetics, Shinilooks. Gelfand, A. E. and Smith, A. F. (1990). Sampling-based approaches to calculating marginal densities, Journal of the American statistical association, 85, 398 409. Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2013). Bayesian Data Analysis (3rd). Chapman and Hall/CRC. Gelman, A. and Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences, Statistical Science, 7, 457 472. Ghahramani, Z. and Beal, M. J. (2001). Propagation algorithms for variational Bayesian learning. In Advances in neural information processing systems, pp. 507 513. Godfrey, K. (1983). Compartment models and their application. Academic Press Inc. London/New York. Ha rdle, W. and Simar, L. (2012). Applied Multivariate Statistical Analysis. Springer. Hastings, W. K. (1970). Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57, 97 109. Hoffman, M. D. and Gelman, A. (2014). The no-u-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15, 1593 1623. Jordan, M. I., Ghahramani, Z., Jaakkola, T. S., and Saul, L. K. (1999). An introduction to variational methods for graphical models. Machine learning, 37, 183 233. Kucukelbir, A., Ranganath, R., Gelman, A., and Blei, D. (2015). Automatic differentiation variational inference in Stan. In Advances in neural information processing systems, pp. 568 576. Kucukelbir, A., Tran, D., Ranganath, R., Gelman, A., and Blei, D. (2017). Automatic differentiation variational inference. Journal of Machine Learning Research, 18, 430 474. Matcuura, K. (2016). Bayesian statistical modeling using Stan and R, Kyoritsu Shuppan Co., Ltd 343. Murphy, K. P. (2012). Machine Learning: A Probabilistic Perspective. MIT Press. Olive, D. J. (2014). Statistical Theory and Inference. Springer. Robert, C. and Casella, G. (2013). Monte Carlo Statistical Methods. Springer.
22 Sun Park, Seongil Jo, Woojoo Lee Shim, C. G. (1994). Pharmacokinetics. Seoul National University Press. Stan Development Team (2017). Stan Modeling Language User s Guide and Reference Manual. Upton, R. A., Thiercelin, J. F., Guentert, T. W., Wallace, S. M., Powell, J. R., Sansom, L., and Riegelman, S. (1982). Intraindividual variability in theophylline pharmacokinetics: statistical verification in 39 of 60 healthy young adults. Journal of Pharmacokinetics and Biopharmaceutices, 10, 123 134. Wainwright, M. J. and Jordan, M. I. (2008). Graphical models, exponential families, and variational inference. Foundations of Trends in Machine Learning, 1, 1 305. Wakefield, J. (2013). Bayesian and Frequentist Regression Methods. Springer. Received October 5, 2020; Revised November 14, 2020; Accepted December 7, 2020
23 약물동태학 모형에 대한 변분 베이즈 방법 박선a, 조성일 1,b, 이우주c a 전북대학교 통계학과, b 인하대학교 통계학과, c 서울대학교 보건대학원 요 약 본 논문에서는 평균장 방법(mean-field methods)을 기반으로 사후 분포(posterior distribution)를 근사하 는 방법인 변분 베이즈 방법(variational Bayes methods)에 대해 소개한다. 특히, 모수들을 실수공간으로 변환 후의 결합 사후분포를 가우시안 분포(Gaussian distribution)들의 곱(product)으로 근사하는 방법인 자동 미분 변분 추론(automatic differentiation variational inference)방법에 대해 자세히 소개하고, 환자에게 약물을 투 여한 후 시간에 따라 약물의 흐름을 파악하는 연구인 약물동태학 모형(pharmacokinetic models)에 적용한다. 소개된 변분 베이즈 방법을 이용하여 자료분석을 실시하고 마코프 체인 몬테 카를로(Markov chain Monte Carlo)방법을 기초로한 자료분석의 결과와 비교한다. 알고리즘의 구현은 Stan을 이용한다. 주요용어: 변분 베이즈, 약물동태학, 자동 미분 변분 추론, 마코프 체인 몬테 카를로, Stan 이 논문은 박선학생 석사논문의 일부를 발췌하였음. 조성일의 연구는 2020년도 정부(과학기술정보통신부)의 재원으로 한국연구재단의 지원을 받아 수행된 기초연구사업임 (NRF-2020R1C1C1A01013338). 교신저자: (22212) 인천광역시 미추홀구 인하로 100, 인하대학교 통계학과. E-mail: statjs@inha.ac.kr