Journal of the Korean Data & Information Science Society 2017, 28(1), 153 161 http://dx.doi.org/10.7465/jkdi.2017.28.1.153 한국데이터정보과학회지 베이지안음이항분기과정을이용한한국메르스발생연구 박유하 1 최일수 2 12 전남대학교통계학과 접수 2016 년 12 월 25 일, 수정 2017 년 1 월 16 일, 게재확정 2017 년 1 월 19 일 요약 전염병확산에대한확률과정모형으로활용되는분기과정은실제데이터를통해모수를추정할수있다는장점이있다. 음이항분포를분기과정의생산분포모형으로적용할수있는데음이항분포를적용하기위해서는평균과산포모수를추정하여야한다. 기존의생물학연구와역학연구분야에서는이를최대우도법을이용하여추정하고있다. 그러나대부분의역학자료의특성상분기과정에서이용되는음이항분포는소표본이어서최대우도추정량의정도를충족시킬수없다. 본논문에서는소표본자료에서좋은통계량의성질을만족한다고알려져있는베이지안을이용하여모수를추정하는방법을제안한다. 2015 년국내메르스사례에베이지안방법을적용하여모수를추정하고사후분포를적합하였다. 그결과어떠한사전분포를가정하더라도안정적으로모수를추정하는것을알수있었다. 추정된산포모수를이용하여분기과정에서의전염병소멸확률을유도하였다. 주요용어 : 베이지안추론, 분기과정, 산포모수, 수학적모델링, 음이항분포, 전염병확산. 1. 연구의배경 2015 년한국에서발생한메르스 (middle east respiratory syndrome coronavirus; MERS-CoV) 는 중동지역에서집중적으로발생한바이러스로사스 (severe acute respiratory syndrome coronavirus; SARS-CoV) 와유사하게주로고열, 기침, 호흡곤란등호흡기증상을일으키는질병이다 (Hwang 과 Oh, 2016). 메르스의전파경로는아직정확히밝혀지지않았지만, 메르스는동물에의한사람으로의 직 간접적인전파또는사람간접촉에의한전파로전염된다. 사우디아라비아에서발생한메르스는낙 타가매개체인것으로추정되며 2015 년한국에서발생한메르스는주로병원내에서메르스감염환자 와의접촉에의해전파되었다. 공기중으로감염되지않고기침이나콧물등호흡기분비물에의해바이 러스가전파되는메르스는특정감염자로부터대다수의 2 차감염자 (offspring; secondary cases; 줄여 서 2 차감염자 ) 가발생하는주요전파사건 (Superspreading Event, SSE) 이다. 메르스는단기간에급속 도로유행하였고한국은전국적으로메르스의공포에휩싸였다. 총 186 명이메르스를확진받았으며이 중다섯명이슈퍼전파자 (superspreader) 이었다 (Korea Centers for Disease Control and Prevention, 2015). 이러한질병확산과관련한연구는생물학, 역학, 질병관리등다양한분야에서관심을가지고연 구되어지고있다. Lloyd-Smith 등, (2005) 은역학모형으로편미분을이용한수학적모형이아닌확 률과정을이용한통계적모형을제시하였다. 주요전파사건으로분류되는여러질병에대한통계적역 본연구는농촌진흥청공동연구사업 ( 과제번호 : PJ01156302) 의지원에의해이루어진것임. 1 (61186) 광주광역시북구용봉로 77, 전남대학교통계학과, 석사과정. 2 교신저자 : (61186) 광주광역시북구용봉로 77, 전남대학교통계학과, 교수. E-mail: ichoi@jnu.ac.kr
154 Yuha Park Ilsu Choi 학모형으로 single-type 분기과정 (branching process) 을적합하였고, 질병의소멸확률 (extinction probability) 을계산하였다. 또한 Alexander (2010) 는돌연변이가고려된병원균의확산모형을 multitype 분기과정을이용하여적합하였고, 이를사용하여소멸확률을계산하였다. 이모형에서분기과정생산분포의모수는지정된값을지닌다. 기초감염재생산수 (basic reproductive number) R 0 는질병에걸릴가능성이있는집단에서감염자에의한 2차감염자수의평균을의미하는데이를보고질병이확산되는정도에대하여대략적으로알수있다. 따라서 R 0 를통해서주요전파사건여부를판단할수있다. 기초감염재생산수를추정하기위한모형은전염병을시간종속적인평균의개념을갖는동역학으로처리한 SIR (susceptible-infetiousrecovery) 모형과확률적추계모형과분기과정을이용한모형으로분류될수있다 (Choi와 Rhee, 2015; Ryu와 Choi, 2015). 질병발생이끝날확률을말하는소멸확률은 R 0 와 k로부터계산된다. Lloyd-Smith 등, (2005) 은분기과정을사용하여 2차감염자의수 Z가 k명일때의확률분포를의미하는생산분포 P (Z = k) 에대하여포아송분포, 기하분포, 음이항분포를가정하고, 주요전파사건에대하여각생산분포에따른소멸확률을구하였다. 소멸확률은질병에대한위험 (risk) 을말해주는지표로소멸확률이높을수록질병이불규칙적으로발생하고일시적인유행을갖는다. 반면낮은소멸확률은질병이만성적으로생존함을의미한다. 여기서환자의감염력 (individual s infectious history) 을설명하는분포로음이항분포는기하분포보다적용범위가넓다고알려져있다 (Burnham과 Anderson, 2002, Lloyd-Smith 등, 2005). 그러므로음이항분포의모수에대한추정량의정도 (precision) 를충족시켜역학자료의모형적합성을높이는것이중요하다고할수있다. 음이항분포의모수는평균 R 0 와산포모수 k이다. Lloyd-Smith 등 (2005) 은음이항분포의산포모수 k를추정하기위해서최대우도법 (maximum-likelihood method) 을이용하였다. 최대우도추정치는소표본오차 (small-sample bias) 가적고효율성의성질이더좋다고알려졌다 (Anscombe, 1950; Pieters 등, 1977). 산포모수 k를추정하는연구는붓스트랩 (bootstrap) 등의방법으로활발히진행되고있지만수치해석적으로구하기때문에소표본추정의편의는모수의참값을과대추정하는단점이있다 (Saha과 Paul, 2005; Pieters, 1977; Piegorsch, 1990). 베이지안방법 (Bayesian method) 은현재주어진관측치에의존하여모수를추정하므로기존방법에비해통계적으로더좋은특성을갖는모수를추정할수있다. 또한베이지안방법은데이터가적을때특히그유용성이좋다고알려져있다 (Hoff, 2009). 본논문에서베이지안방법을사용하여음이항분포를생산분포로하는분기과정의모수를추정하는방법을제시한다. 또한 2015년한국에서발생한메르스자료에이를적합시켜소멸확률을구하였다. 2015년한국의메르스를비롯하여대부분의역학자료는소표본이므로분기과정을이용한모형을적합할시에기존의방법으로적합시킨모형보다베이지안방법으로적합시킨모형이더타당하다고할수있다. 2. 연구자료본연구에서는 2015년한국메르스의분기과정을파악하기위해서보건복지부질병관리본부와 WHO에서조사되어공개된자료를활용하였다. 2015년 5월 20일확진된한국의첫메르스환자는증상발현후약 10일동안여러병 의원을다니며가족과다른환자, 의료진과접촉하여다수의 2차환자를유발하였다. 한달이채되지않는기간동안한국에서메르스확진을받은전체환자의수는 186명이다. 본연구에서는 2015년 5월초부터 6월말까지의두달동안감염된총 186명의환자를대상으로메르스감염경로를추적하였다. 감염경로가확실하지않은 17명의환자를제외한 169명의환자에대한자
A study on MERS-CoV outbreak in Korea using Bayesian negative binomial branching processes 155 료로부터 2015 년 5 월 20 일발생한한국의첫환자로부터의감염경로를파악하여분기과정을구성하면 Figure 2.1 과같다. Figure 2.1 Branching process of MERS in Korea (2015) 3.1. 베이지안음이항분기과정모형 3. 베이지안음이항분기과정모형 분기과정은 Galton-Watson 분기과정이라고도하며 1873 년에 Francis Galton 에의해서처음으로제 안된것이다 (Kendall, 1966). 본논문에서는 single-type discrete time 분기과정을 2015 년한국메르 스자료에적용하였다. single-type discrete time 분기과정은같은세대 (generation) 에서질병을전파 할가능성이있는모든환자는다음세대에서동일한성질의질병을갖는 2 차감염자를생산하고이때, 2 차감염자의수 Z (Z = 0, 1, 2,, n) 는서로독립이고동일한생산분포를따른다고가정한다. 생산 분포는 2 차감염자의수 Z 가 k 명일때의확률분포를의미하며, P (Z = k) = p k 로주어진다. i 번째 개체의 2 차감염자의수를 Z i 이라고하고, n 번째세대에서감염자수의총합을 X n 라고하면, X n 은 1, 2,, Z n 1 감염자에게서전염된 2 차감염자수의총합이다 (Kimmel, 2002; Haccou, 2005; Pinsky, 2011). 즉, X n = X n 1 i=1 Z i. (3.1) 감염재생산수 ν 는특정한환자에의해서감염된 2 차감염자의감염기간을의미하며환자의감염력 에서의모든변동을나타낸다 (Lloyd-Smith 등, 2005). 또한감염재생산수 ν 의값은모평균이기초감 염재생산수 R 0 인연속확률분포로부터추정할수있다. 주요전파사건 (superspreading events) 은감염 재생산수 ν 의분포가오른쪽꼬리분포의형태로나타난다는특징을갖는다 (Dye 등, 2003; Anderson, 2004; Lloyd-Smith 등, 2005).
156 Yuha Park Ilsu Choi 본논문에서는생산분포를음이항분포로가정한다. 우선, 전염병의확률적효과가포아송과정을 따른다고가정하면 2 차감염자의수 Z 는평균이감염재생산수 ν 인포아송혼합분포 (poisson mixture distribution) Z P oisson(v) 를따른다고할수있다 (Lloyd-Smith 등, 2005). 여기서, 감염재생산수 ν 의분포를감염재생산수 ν 가평균이 R 0 이고산포모수 (dispersion parameter) 가 k 인감마분포를따를때즉, v Gamma (k, R 0/k) 라고한다면생산분포는결국음이항분포 ( Z NB ( (1 + R0/k) 1, k )) 로주어진다 (Patil, 1970; Pielou, 1977; Pinsky 과 Karlin, 2011). 분기과정에서생산분포는확률생성함수 (probability generating function, pgf) 로나타낸다. 그러므 로베이지안음이항분기과정모형은다음과같은음이항분포를생산분포로갖는분기과정으로표현된 다. 즉, g(s) = E(s Zn ) = ( q 1 ps ) n, s 1. (3.2) 이때, p = (1 + R 0/k) 1 이며, g 0(s) = s, g 1(s) = g(s),, g n+1(s) = g(g n(s)), n = 1, 2, 3, 을 만족한다. 기초감염재생산수 R 0 는 Z 의평균이며 g (1) 로나타낼수있기때문에모수추정에있어기초감염재 생산수 R 0 는표본의평균으로간단히추정할수있다. 그러나산포모수 k 의경우간단히추정하기어 려우므로본연구에서는베이지안추론을이용하였다. 베이즈정리 (Bayes rule) 에의해서실제관측된관측치 Z = (z 1, z 2,, z n) 는평균이 v 인포아송 분포를따르고 ν 는평균이 R 0 이고산포모수가 k 인감마분포를따른다고하자. 이때사후분포는 Z ν P oisson(ν), (3.3) ν Gamma (k, R 0/k). (3.4) f(k Z) f (Z ν) f (ν k, R 0/k) f (k, R 0/k) (3.5) n ν z i e ν R 0/k z i! Γ(k) (R0ν/k)k 1 e R 0ν/k f (k, R 0/k). (3.6) i=1 그러므로 p = (1 + R 0/k) 1 일때, f(z ν)f (ν k, R 0/k) 로나타낼수있고최종적으로식 (3.5) 는 가된다. f(k Z) ( ) n k + zi 1 (1 p) k (p) z i (3.7) k i=1 ( ) n k + zi 1 (1 p) k f (k, R 0/k) (3.8) k i=1 위에서알수있는것처럼모수에대한베이지안추론은모수의사후분포에근거한다. 식 (3.5) 에서 모수의사후분포는자료가주어졌을때모수의사전분포에의존한다는것을알수있다. 베이지안방법 에서사전분포는자료를얻기전에모수에대한정보를나타내기때문에적절한사전분포의선택이중 요하다. 사전분포는두가지의경우즉, 과거의경험이나이해에따른주관적사전정보를요약한주
A study on MERS-CoV outbreak in Korea using Bayesian negative binomial branching processes 157 관적사전분포 (subjective prior) 와사전정보가확실하지않거나없을때사용하는객관적사전분포 (objective prior) 또는무정보적사전분포 (non-informative prior) 로구분된다. 각각의경우에서연구자의선택에따라매우다양한분포를가정할수있다 (Hogg, 2009). 본연구에서는실제자료에근거하여, 모수 k에대해최소한의사전정보만을이용한무정보적사전분포로일양분포 (uniform prior) 를고려하였다. 일양사전분포의모수는 k의범위에기초하여 (0, 20) 으로가정하였다. 또한 Lloyd-Smith 등 (2005) 의연구결과에기초하여, 주관적사전분포로지수분포 (exponential prior) 와로그정규분포 (log-normal prior) 도고려하였다. 2015년한국메르스자료의기본감염재생성수 R 0 는 1.98이므로지수사전분포의모수를 2로가정하였으며로그-척도모수는 0이고형상모수는 1인로그정규사전분포를가정하였다. 3.2. 시뮬레이션과모형진단사후밀도함수를계산하는과정에서필요한복잡한형태의함수적분을하기위해서베이지안추론에서대표적으로사용하는시뮬레이션방법인마코브체인몬테칼로 (Markov Chain Monte Carlo; MCMC) 알고리즘을이용할수있다. 본연구에서는 MCMC 기법의일종인깁스표집기 (Gibbs sampler) 를사용하였다. 깁스표집기를위한 v의사후분포는앞절에서설명한것처럼식 (3.8) 과같고이에대한모수 k의분포는일양사전분포, 지수사전분포, 로그정규사전분포를가정하였다. 각사전분포별로깁스표집기에서 10,000의반복샘플링을하였으며 500개를번인 (burn-in) 하였다. 본논문에서는난수를생성하여모수를추정한후, 수렴이잘되었는지를판단하기위하여 Gelman- Rubin 통계량을사용하였다. 체인 (Chain) 은 10개로설정하였고각각의체인에대한샘플링을 10,000번하였다. Gelman-Rubin 통계량은 1과가까운값을가질수록잘수렴한다 (Gelman 과 Rubin, 1992; Monahan, 2001; Jo와 Lee, 2016). 3.3. 소멸확률본연구에서는 2015년한국메르스에대해서질병의위험정도를분석하기위해소멸확률을구하였다. X n 의확률생성함수가 g n(s) 라면, n번째세대에서질병의소멸확률 (extinction probability) P (X n = 0) 은 g n(0) 으로나타낼수있다. 여기서이론적으로 g (1) < 1이면, 질병은소멸된다. 만약 g (1) = 1이면, 질병은소멸되지만소멸확률이정해지지않게된다. 또한만약 g (1) > 1이면, 질병은확산된다. 따라서 g (1) = R 0 1일때, 생산분포의확률생성함수를 g(s) 라고하면소멸확률을 g(s) = s로정의한다. 이러한소멸확률은수학적으로구하기어렵기때문에일반적으로수치해석적인방법을사용하여유도한다. 4. 자료분석결과 2015년한국메르스자료에대하여산포모수 k를베이지안추론을이용하여추정한결과는 Table 4.1이다. Figure 4.1은시뮬레이션을통해생성된표본들의수렴성을나타내는표본의경로그림이다. 고려된사전분포들이모두수렴함을알수있다. 그러나로그정규분포일때의진폭이다른두사전분포에비해넓어수렴에는불리한결과를갖고있음을알수있다. k exp(2) 일경우, k의중앙값을 1.409, k의표준편차를 0.093으로추정하여가장나은성능을보였다. 반면 k lognormal (0, 1) 일경우, k의중앙값을 1.414로 k의표준편차를 0.094로추정하였다. 그러나세가지분포를가정한경우모두추정된 k의중앙값은비슷하였고 k의표준편차의추정값또한별다른차이가없었다.
158 Yuha Park Ilsu Choi Table 4.1 Bayesian estimates based on prior distributions prior mean sd MC error 2.5% median 97.5% unif(0, 20) 1.430 0.093 0.0018 1.252 1.428 1.616 exp(2) 1.411 0.093 0.0021 1.237 1.409 1.596 lognormal(0, 1) 1.417 0.094 0.0018 1.242 1.414 1.608 Figure 4.1 MCMC trace plot according to uniform, exponential and log-normal prior Figure 4.2는체인 (chain) 을 10개로설정하여구한 Gelman-Rubin 통계량을나타낸그림과통계량이다. 푸른색은체인내변이 (within-chain variance) 의추정치를, 초록색은체인간변이 (betweenchain variance) 의추정치를나타내며붉은색은체인내변이의추정치와체인간변이의추정치의비 (ratio) 를나타낸다. 이러한비를 Gelman-Rubin 통계량이라고하며자료분석결과 1로수렴하는것을알수있다. 따라서어떠한사전분포를가정하더라도안정적으로모수가수렴하는것을알수있다. Figure 4.2 Gelman-Rubin plot Figure 4.3은일양분포, 지수분포와로그정규분포사전확률분포에따른베이지안전염병소멸확률이다. 확률생성함수 y = g(s) 와직선 y = s가교차되는것을구하는것으로사전확률분포에따른베이지안전염병소멸확률을유도하면 Table 4.2와같다. 추정된산포모수 k의중앙값을가진세가지의사전확률분포에서뉴턴-랩슨알고리즘을사용하여각각 5,000개의표본을추출하여소멸확률을구하였
A study on MERS-CoV outbreak in Korea using Bayesian negative binomial branching processes 159 다. 소멸확률은일양사전분포를가정하는경우 0.445 로가장크게나타났고지수사전분포를가정하 는경우 0.434 로가장낮게나타났다. Table 4.2 Bayesian extinction probability based on prior distributions prior estimate of k extinction probability unif(0, 20) 1.428 0.445 exp(2) 1.409 0.434 lognormal(0, 1) 1.414 0.439 Figure 4.3 Bayesian extinction probability according to uniform, exponential and log-normal prior 5. 결론본연구에서는분기과정에서생산분포가음이항분포를따를때, 모수의재구성 (reparameterization) 을통하여베이지안추론을유도하였다. 산포모수 k의사후분포를유도하기위해서첫째, 일양사전분포를따르는경우, 둘째, 지수사전분포를따르는경우, 셋째, 로그정규사전분포를따르는경우를가정하였다. 본논문에서사용된베이지안음이항분기과정은세가지의장점을가진다. 첫째, 2015년한국메르스의실제자료를이용하여모형에적용한결과, 베이지안방법을이용하면산포모수 k의추정을안정적으로할수있음을알수있었다. 세가지의사전분포에서산포모수 k의베이지안추정값의차이가많지않았으며, 추정값의표준편차에서도많은차이가나지않았다. 또한 Gelman-Rubin 통계량을사용하여모형진단을한결과통계량값이 1에가까워모수의수렴성이확인되었다. 둘째, 베이지안추론을이용하면기존추정치의문제점으로지적되는모수의참값을과대추정하는단점을해결하고수치해석적으로구하는최대우도법의이론적오류를해결할수있을것으로판단된다. 셋째, 분기과정을사용하면질병의소멸확률을계산할수있다는장점이있다. 본연구에서는베이지안음이항분기과정을이용해서 2015년한국메르스자료의위험정도를판단하기위하여소멸확률을계산하였다. 세가지의사전분포가정에서구한소멸확률은 43% 에서 44% 로, 2015년의한국메르스는낮은소멸확률을갖는특이한현상을나타내었다. 따라서한국메르스는숙주인낙타가없음에도불구하고전염력이항상존재하는질병으로발전할수있는가능성을가지고있음을알수있다. 그러나본논문에서는기존의최대우도법추정치와비교를실시하지않았다. 이후의연구에서는이와관련된연구를진행할예정이며주요전파사건과역학적관련성이있는질병에대한다양한자료를수집하여본연구에서사용한베이지안음이항분기과정을역학연구에일반적으로활용할수있음을연구할예정이다.
160 Yuha Park Ilsu Choi References Alexander, H. K. (2010). Modelling pathogen evolution with branching processes, Master Thesis, Queen s University, Ontario. Anderson, R. M., Fraser, C., Ghani, A. C., Donnelly, C. A., Riley, S., Ferguson, N. M., Leung, G. M., Lam, T. H. and Hedley, A. J. (2004). Epidemiology, transmission dynamics and control of SARS: The 2002-2003 epidemic. Philosophical Transactions of the Royal Society B, 359, 1091-1105. Anscombe, F. J. (1950). Sampling theory of the negative binomial and logarithmic series distributions. Biometrika, 37, 358-382. Choi, I. and Rhee, S-S. (2015). A transmission distribution estimation for real time Ebola virus disease epidemic model. Journal of the Korean Data & Information Science Society, 26, 161-168. Dye, C. and Gay, N. (2003). Modeling the SARS epidemic. Science, 300, 1884-1885. Gelman, A. and Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences (with discussion). Statistical Science, 7, 457-511. Haccou, P., Jagers, P. and Vatutin, V. A. (2005). Branching processes: Variation, growth, and extinction of populations, Cambridge University Press, New York. Hoff, P. D. (2009). A first course in Bayesian statistical methods, Springer Dordrecht Heidelberg London, New York. Hwang, S. and Oh, C. (2016). Estimation of the case fatality ratio of MERS epidemics using information on patients severity condetion. Journal of the Korean Data & Information Science Society, 27, 599-607. Jo, S., and Lee, J. (2016). Density estimation of summer extreme temperature over South Korea using mixtures of conditional autoregressive species sampling model. Journal of the Korean Data & Information Science Society, 27, 1155-1168. Kendall, D. G. (1966). Branching processes since 1873. Journal of London Mathematics Society, 41, 386. Kimmel, M. and Axelrod, D. E. (2002). Branching processes in biology, Springer-Verlag Inc., New York. Korea Centers for Disease Control and Prevention. (2015). Middle east respiratory syndrome coronavirus outbreak in the Republic of Korea. Osong Public Health and Research Perspectives, 6, 269-278. Lloyd-Smith, J. O., Schreiber, S. J., Kopp, P. E., Getz, W. M. (2005). Superspreading and the impact of individual variation on disease emergence. NATURE, 438, 355-359. Monahan, J. F. (2001). Numerical methods of statistics, Cambridge University Press, New York. Patil, G. P. (1970). Random counts in models and structures, Pensylvania State University Press, University Park PA. Piegorsch, W. W. (1990). Maximum-likelihood estimation for the negative binomial dispersion parameter. Biometrics, 46, 863-867. Pielou, E. C. (1977). Mathematical Ecology, Wiley, New York. Pieters, E. P., Gates, C. E., Matis, J. H. and Sterling, W. L. (1977). Small sample comparison of different estimators of negative binomial parameters. Biometrics, 33, 718-723. Pinsky, M. A. and Karlin, S. (2011). An Introduction to stochastic modeling, 4th Ed., Academic Press, San Diego. Ryu, S. and Choi, B. (2015). Development of epidemic model using the stochastic method. Journal of the Korean Data & Information Science Society, 26, 301-312. Saha, K. and Paul, S. (2005). Bias-corrected maximum likelihood estimator of the negative binomial dispersion parameter. Biometrics, 61, 179-185.
Journal of the Korean Data & Information Science Society 2017, 28(1), 153 161 http://dx.doi.org/10.7465/jkdi.2017.28.1.153 한국데이터정보과학회지 A study on MERS-CoV outbreak in Korea using Bayesian negative binomial branching processes Yuha Park 1 Ilsu Choi 2 12 Department of Statistics, Chonnam University Received 25 December 2016, revised 16 January 2017, accepted 19 January 2017 Abstract Branching processes which is used for epidemic dispersion as stochastic process model have advantages to estimate parameters by real data. We have to estimate both mean and dispersion parameter in order to use the negative binomial distribution as an offspring distribution on branching processes. In existing studies on biology and epidemiology, it is estimated using maximum-likelihood methods. However, for most of epidemic data, it is hard to get the best precision of maximum-likelihood estimator. We suggest a Bayesian inference that have good properties of statistics for small-sample. After estimating dispersion parameter we modelled the posterior distribution for 2015 Korea MERS cases. As the result, we found that the estimated dispersion parameter is relatively stable no matter how we assume prior distribution. We also computed extinction probabilities on branching processes using estimated dispersion parameters. Keywords: Bayesian inference, branching processes, dispersion parameter, epidemic emergence, mathematical modeling, negative binomial distribution. This study was carried out with the support of the Research Program of Rural Development Administration, Republic of Korea (Project No. PJ01156302). 1 Master degree, Department of Statistics, Chonnam National University, Gwangju 61186, Korea. 2 Corresponding author: Professor, Department of Statistics, Chonnam National University, Gwangju 61186, Korea. E-mail: ichoi@jnu.ac.kr