處理零計數(shù)過多數(shù)據(jù)的兩部模型方法介紹*

哈爾濱醫(yī)科大學(xué)衛(wèi)生統(tǒng)計學(xué)教研室（150081）

姜 博 門志紅 劉匆提 劉 艷△

·綜述·

處理零計數(shù)過多數(shù)據(jù)的兩部模型方法介紹*

哈爾濱醫(yī)科大學(xué)衛(wèi)生統(tǒng)計學(xué)教研室（150081）

姜 博 門志紅 劉匆提 劉 艷△

在公共衛(wèi)生、社會科學(xué)等領(lǐng)域的調(diào)查研究中，定量數(shù)據(jù)可能會包含過多零值，定性數(shù)據(jù)也有可能存在著零膨脹［1］，如果仍按照一般模型（如泊松回歸或線性回歸）的數(shù)據(jù)分布假設(shè)進(jìn)行分析，會導(dǎo)致參數(shù)估計產(chǎn)生偏倚［1-2］。針對此類問題，兩部模型法是一種比較常用的擬合模型方式。目前，兩部模型法包括應(yīng)用于定性數(shù)據(jù)的 Hurdle模型、零膨脹模型（zero-inflated model）等，以及應(yīng)用于定量數(shù)據(jù)的半連續(xù)數(shù)據(jù)（semicontiunuous data）兩部模型等［3］。

定性數(shù)據(jù)兩部模型

對于定性數(shù)據(jù)的零膨脹問題，可以采用混合泊松回歸模型、birth process模型、兩部模型等［4］。兩部模型的基本思想是將數(shù)據(jù)分成兩個部分，第一部分考慮事件是否發(fā)生，第二部分考慮事件發(fā)生的次數(shù)，主要包括 hurdle模型［5］與零膨脹模型［6］。如 Neelon等調(diào)查某一地區(qū)內(nèi)發(fā)生急診科就醫(yī)的情況，因其中大部分被調(diào)查者急診科就醫(yī)為零計數(shù)，故使用hurdle模型進(jìn)行分析［7］；Rose等使用零膨脹模型分析疫苗不良事件發(fā)生情況［8］。兩種模型主要的區(qū)別在于零值來源的假設(shè)處理不同［9］，Hurdle模型更適合處理數(shù)據(jù)的零計數(shù)僅來自于抽樣的零，而零膨脹模型更適合處理數(shù)據(jù)的零計數(shù)來自結(jié)構(gòu)的零與抽樣的零，實際應(yīng)用中還要根據(jù)研究目的、專業(yè)解釋及數(shù)據(jù)的具體分布來選擇合適的兩部模型［10］。

1.Hurdle模型

Mullahy對Hurdle模型進(jìn)行了深入討論［5］，第一部分考慮零計數(shù)是否發(fā)生，服從二項分布假設(shè)，連接函數(shù)一般為 logit、clog、probit等［9］；第二部分再對非零計數(shù)部分進(jìn)行分析，一般假設(shè)其服從泊松分布，擬合泊松回歸模型，即泊松hurdle模型。此外，第二部分還有其他的分布假設(shè)類型，如負(fù)二項分布［10］、廣義泊松分布［11］等。Hurdle模型參數(shù)估計是將兩個部分作為相互獨立的模型進(jìn)行估計，其具體方法包括多種，較常用的是極大似然法，也有使用廣義估計方程（generalized estimating equations，GEE）方法分別對兩部分進(jìn)行參數(shù)估計［12］。

以泊松hurdle模型為例，其第一部分為p（y＝0｜x）＝π，0≤π≤1，第二部分為p（y＞0｜x）＝（1-中π為發(fā)生零計數(shù)的概率，μ為事件發(fā)生的平均次數(shù)，y為事件發(fā)生次數(shù)［5］；其參數(shù)估計的對數(shù)似然函數(shù)為：

2.零膨脹模型（zero-inflated model）

Lambert［6，10］首次提出零膨脹泊松模型，模型的假設(shè)是數(shù)據(jù)服從一種混合分布［3］，數(shù)據(jù)分為兩個部分，同時數(shù)據(jù)中的“0”也被分成兩種組成部分，第一部分仍然考慮是否發(fā)生零計數(shù)情況，假設(shè)服從二項分布，連接函數(shù)一般為logit，該部分的零值稱為不可能發(fā)生的零或結(jié)構(gòu)的零，即不存在發(fā)生可能所導(dǎo)致的零值；第二部分考慮可能發(fā)生事件的情況，即抽樣的零，也就是有可能發(fā)生但并未發(fā)生事件，常用的模型假設(shè)為泊松分布，負(fù)二項分布［10］、廣義泊松分布［13］、廣義冪級數(shù)分布（generalized power series distribution）［14］等。零膨脹模型的參數(shù)估計依然較常采用極大似然法，可使用最大期望算法（expectation maximization algorithm，EM）［6］、牛頓迭代法（Newton-raphson algorithm）［15］進(jìn)行極大似然運算，同樣也可以使用廣義估計方程方法估計相應(yīng)參數(shù)［16］，或應(yīng)用貝葉斯方法進(jìn)行參數(shù)估計［17］。

以零膨脹泊松模型［6，10］為例，其第一部分為 p意義與 hurdle模型中相同；其對數(shù)似然函數(shù)表示為［6，10］exp（x1β2）］。

3.模型發(fā)展

伴隨著定性數(shù)據(jù)兩部模型方法的不斷發(fā)展，兩部模型的應(yīng)用也不斷擴展。已有一些研究在零膨脹模型的基礎(chǔ)上提出了半?yún)?shù)零膨脹模型［18］、修正零膨脹泊松模型［1］等模型，還有一些研究也提出適合于其他類型數(shù)據(jù)應(yīng)用的模型，如縱向定性數(shù)據(jù)的零膨脹增長曲線模型（zero-inflated growth curve model）［19］、層次結(jié)構(gòu)數(shù)據(jù)的多水平零膨脹泊松模型（multi-level zero-inflated poisson regression modelling）［20］、空間泊松 hurdle模型（spatial poisson hurdle model）［7］、時空 hurdle模型（spatiotemporal hurdle models）［21］等等。

定量數(shù)據(jù)兩部模型

在大型隊列研究或橫斷面調(diào)查中，醫(yī)療費用經(jīng)常呈現(xiàn)為一種半連續(xù)數(shù)據(jù)，即在調(diào)查人群中存在大量的調(diào)查對象未發(fā)生醫(yī)療行為即不產(chǎn)生費用，發(fā)生醫(yī)療行為調(diào)查對象的費用呈正偏態(tài)分布。對于這類半連續(xù)數(shù)據(jù)，過多的零值導(dǎo)致費用分布右偏嚴(yán)重，一般傳統(tǒng)模型假設(shè)分布類型難以滿足這類半連續(xù)數(shù)據(jù)［22］，可以選用Tobit模型［23］、樣本選擇（sample selection）模型［24］、Cox比例風(fēng)險回歸模型［25］、兩部模型［26］進(jìn)行數(shù)據(jù)分析。如Finkelstein等使用美國MEPS（medical expenditure panel surveys）數(shù)據(jù)，應(yīng)用兩部模型分析肥胖與醫(yī)療費用的關(guān)系［27］；Bock等使用隊列人群的橫斷面調(diào)查數(shù)據(jù)，應(yīng)用兩部模型分析德國老年人口自費健康醫(yī)療服務(wù)的不公平性問題［28］。

1.模型組成

與定性數(shù)據(jù)的hurdle模型類似，定量數(shù)據(jù)的兩部模型依然是將半連續(xù)數(shù)據(jù)分成兩個部分，第一部分將應(yīng)變量作為二項分布處理，如是否發(fā)生醫(yī)療行為，構(gòu)建二項分布概率模型，常用logistic、probit回歸模型分析是否發(fā)生醫(yī)療行為的影響因素；第二部分對發(fā)生醫(yī)療行為即大于零的數(shù)據(jù)部分，一般常假設(shè)費用滿足對數(shù)正態(tài)分布［29］，進(jìn)行數(shù)據(jù)轉(zhuǎn)換后擬合廣義線性回歸模型，分析醫(yī)療費用的影響因素，gama分布也可以作為第二部分模型的分布假設(shè)［30］。當(dāng)假設(shè)費用數(shù)據(jù)服從對數(shù)正態(tài)分布時，如果存在個體真實費用為零的情況時，由于零值不能進(jìn)行對數(shù)轉(zhuǎn)化，若將其歸為第二部分，則軟件運算時會將其當(dāng)作缺失值而忽略。

2.參數(shù)估計

半連續(xù)數(shù)據(jù)兩部模型的參數(shù)估計方法較多，選擇也較為靈活，主要是根據(jù)調(diào)查目的及各部分所選用的模型而定，較常用的依然是極大似然法［3］。但隨著兩部模型的應(yīng)用拓展，不同種類兩部模型不斷被開發(fā)出來，其參數(shù)估計方法也各不相同，如MK Olsen與JL Schafer應(yīng)用兩部模型分析縱向數(shù)據(jù)，擬合兩部隨機效應(yīng)模型，使用極大似然法估計固定效應(yīng)，基于高階拉普拉斯法的類Fisher評分法進(jìn)行方差估計［31］；Neelon等使用貝葉斯兩部潛類別模型（two-part latent class model），應(yīng)用馬氏鏈蒙特卡羅（MCMC）法擬合兩部模型［32］。

以第一部分使用probit回歸模型、第二部分使用對數(shù)變換后線性模型為例，擬合兩步法模型，其各部函數(shù)及似然函數(shù)表示為［29］：第一部分probit回歸模型為yi＝β1xi+e1i，e1i～N（0，1），第二部分對數(shù)轉(zhuǎn)化線性模型為 log（yi｜yi＞0）＝β2xi+e2i，e2i～N（0，σ2），表示第i個觀測。

3.模型預(yù)測

半連續(xù)數(shù)據(jù)兩部模型可以用來預(yù)測個體的醫(yī)療費用情況，具體方法是通過第一部分的模型預(yù)測個體可能發(fā)生醫(yī)療行為的概率，再通過第二部分模型預(yù)測所產(chǎn)生費用的期望，兩個部分模型的估計值相乘即可得出個體醫(yī)療費用的估計值，其函數(shù)可以表示為E（yi｜xi）＝prob（yi＞0｜xi）·E（yi｜xi，yi＞0）［33］。通常，第一部分產(chǎn)生費用的概率估計值采用probit回歸或logistic回歸模型進(jìn)行計算，而第二部分具體費用的估計值則可根據(jù)數(shù)據(jù)的分布選擇合適的計算方式。若假設(shè)第二部分?jǐn)?shù)據(jù)服從對數(shù)正態(tài)分布，可依據(jù)最小二乘法（OLS）估計相應(yīng)統(tǒng)計量，再對估計值取指數(shù)，表示為E（yi｜χi，yi＞0）＝exp（β2xi+σ2／2）；若應(yīng)用對數(shù)轉(zhuǎn)化后，誤差項仍無法滿足正態(tài)分布，使用最小二乘法進(jìn)行估計將出現(xiàn)偏倚，故Duan提出一種非參數(shù)的估計方法，即 Smearing估計法［34］。Smearing估計法對數(shù)據(jù)分布沒有特定假設(shè)，僅要求誤差項獨立同分布，模型的估計值可表示為［29］

最小二乘法、Smearing估計法、GLMs模型三種估計方式各有相應(yīng)的適用條件，需根據(jù)數(shù)據(jù)分布的具體情況選擇合適的方法，若數(shù)據(jù)轉(zhuǎn)換后符合相應(yīng)分布（如對數(shù)正態(tài)分布）或峰度較高，則OLS法估計更為準(zhǔn)確；若不符合相應(yīng)分布，可選用Smearing法；轉(zhuǎn)換后數(shù)據(jù)峰度較低時，可使用 GLMs［3］。

4.模型發(fā)展

兩部模型的簡單與靈活性決定其被廣泛應(yīng)用于半連續(xù)數(shù)據(jù)分析［26，36］。1981年 Manning等將兩部模型應(yīng)用于衛(wèi)生經(jīng)濟(jì)領(lǐng)域，研究醫(yī)療保險費用［37］。為滿足多種數(shù)據(jù)分布類型的發(fā)展，兩部模型不斷拓展，第二部分模型已提出廣義 gamma分布、box-cox變換等［38］；處理縱向數(shù)據(jù)、層次結(jié)構(gòu)數(shù)據(jù)的兩部模型已提出兩部隨機效應(yīng)模型［26］、多水平兩部模型［39］、基于貝葉斯方法的雙變量兩部模型 （bivariate two-part model）［40］等；甚至在兩部模型理論的基礎(chǔ)上，提出三部模型［2］、四部模型［29］等等。

相關(guān)軟件

在統(tǒng)計學(xué)軟件方面，R軟件可以使用pscl等程序包［41］，SAS軟件可以通過 GENMOD、NLMIXED等過程，Stata軟件可以通過 ZIP、HPLOGIT等命令、Win-BUGS軟件［42］通過貝葉斯算法擬合構(gòu)建定性數(shù)據(jù)兩部模型；SAS軟件通過GLIMMIX和MCMC（貝葉斯算法）等過程、Stata通過GLM等命令、WinBUGS軟件［43］（貝葉斯算法）構(gòu)建定量數(shù)據(jù)兩部模型。

結(jié) 論

兩部模型可有效地解決公共衛(wèi)生、社會科學(xué)等領(lǐng)域調(diào)查研究中出現(xiàn)的零計數(shù)過多問題，其靈活的數(shù)據(jù)分布類型假設(shè)與參數(shù)估計方式，使其擁有較為廣泛的適用范圍。定性數(shù)據(jù)兩部模型的使用選擇，要充分考慮數(shù)據(jù)中零計數(shù)的來源方式、數(shù)據(jù)分布類型、參數(shù)估計方式、研究目的及專業(yè)知識等方面從而有效減小偏倚；定量數(shù)據(jù)兩部模型使用時，要選擇合適的數(shù)據(jù)分布類型、參數(shù)估計方式及模型預(yù)測的計算方式。目前，仍有較多研究和項目關(guān)注于兩部模型法的計算與使用，使其在零計數(shù)過多問題的處理上繼續(xù)保持較大的應(yīng)用價值及推廣意義。

［1］郭念國.零膨脹泊松模型的改進(jìn)在零次索賠建模中的應(yīng)用.統(tǒng)計與信息論壇，2010，25（7）：22-25.

［2］葉玲瓏.葉玲瓏.基于兩部模型的家庭醫(yī)療需求與消費結(jié)構(gòu)研究.廈門：廈門大學(xué)，2014.

［3］Neelon B.Two-Part Models for Zero-Modified Count and Semicontinuous Data.Duke University，2013.

［4］Ridout M，Demétrio CG，Hinde J.Models for count data with many zeros／／Proceedings of the XIXth International Biometric Conference，1998，19：179-192.

［5］Mullahy J.Specification and testing of some modified count data models.J Econom，1986，33（3）：341-365.

［6］Lambert D.Zero-inflated Poisson regression，with an application to defects in manufacturing.Technometrics，1992，34（1）：1-14.

［7］B Neelon PG，Loebs PF.A spatial Poisson hurdle model for exploring geographic variation in emergency department visits，2013，176（2）：389-413.

［8］CE Rose SWM，Wannemuehler KA.On the use of zero-inflated and hurdle models for modeling vaccine adverse event count data，2006，16（4）：463-481.

［9］曾平，趙晉芳，劉桂芬.居民就診次數(shù)的Hurdle模型分析.中國衛(wèi)生統(tǒng)計，2010，27（6）：603-605.

［10］原靜，劉桂芬，薛玉強.零膨脹計數(shù)資料模型選擇與比較.中國衛(wèi)生統(tǒng)計，2011，28（4）：354-356，360.

［11］SE Saffari RA，Greene W.Investigating the impact of excess zeros on hurdle-generalized Poisson regression model with right censored count data.Statistica Neerlandica，2013，67（1）：67-80.

［12］MJDobbie AHW.Theory＆Methods：Modelling Correlated Zero-inflated Count Data.2001，43（4）：431-444.

［13］H Joe RZ.Generalized Poisson distribution：the property of mixture of Poisson and comparison with negative binomial distribution.Biometrical Journal，2005，47（2）：219-229.

［14］Statistics YZF-Ci.Model Selection of Zero-Inflated Generalized Power Series Distribution with Missing Responses.Theory and Methods，2012，41（6）：1013-1028.

［15］閆鳴宇.零膨脹模型的若干問題研究.蘇州：蘇州大學(xué)，2011.

［16］DB Hall ZZ.Marginal models for zero inflated clustered data.Statistical Modelling，2004，4（3）：161-180.

［17］Statistics JR-Ci.Bayesian analysis of zero-inflated distributions.Theory and Methods，2003，32（2）：281-289.

［18］KF Lam HX，Cheung YB.Semiparametric Analysis of Zero-Inflated Count Data.Biometrics，2006，62（4）：996-1003.

［19］Liu H.Growth curve models for zero-inflated count data：An application to smoking behavior，2007，14（2）：247-279.

［20］AH Lee KW，JA Scott KY.Multi-level zero-inflated Poisson regression modelling of correlated count data with excess zeros，2006，15（1）：47-61.

［21］B Neelon HHC，Ling Q.Spatiotemporal hurdle models for zero-inflated count data：Exploring trends in emergency department visits.Stat Methods Med Res，2014：0962280214527079.

［22］Liu L，Strawderman RL，Cowen ME，et al.A flexible two-part random effects model for correlated medical costs.JHealth Econ，2010，29（1）：110-23.

［23］Tobin J.Estimation of relationships for limited dependent variables.Econometrica：journal of the Econometric Society，1958，26（1）：24-36.

［24］Madden D.Sample selection versus two-part models revisited：The case of female smoking and drinking.Journal of Health Economics，2008，27（2）：300-307.

［25］A Basu WGM，Mullahy J.Comparing alternative models：log vs Cox proportional hazard.Health Econ，2004，13（8）：749-765.

［26］Tian L，Huang J.A two-part model for censored medical cost data.Stat Med，2007，26（23）：4273-92.

［27］EA Finkelstein JGT，JW Cohen WD.Annual medical spending attributable to obesity：payer-and service-specific estimates.Health Aff（Millwood），2009，28（5）：w822-w831.

［28］JO Bock HM，Brenner H.Inequalities in out-of-pocket payments for health care services among elderly Germans-results of a populationbased cross-sectional study.Int JEquity Health，2014，13（1）：3.

［29］N Duan JPN，CN Morris WGM.A comparison of alternative models for the demand for medical care：Rand Corporation，1982：15-32.

［30］周虹，余松林，向蕙云.兩部模型的基本原理與SAS GLIMMIX過程的配合.中國衛(wèi)生經(jīng)濟(jì)，2008，27（10）：19-22.

［31］MK Olsen JLS.A two-part random-effects model for semicontinuous longitudinal data.2001，96（454）：730-745.

［32］Neelon B，O'Malley AJ，Normand SL.A bayesian two-part latent class model for longitudinal medical expenditure data：assessing the impact of mental health and substance abuse parity.Biometrics，2011，67（1）：280-289.

［33］ Buntin MB，Zaslavsky AM.Too much ado about two-part models and transformation？Comparing methods of modeling Medicare expenditures.JHealth Econ，2004，23（3）：525-542.

［34］Duan N.Smearing estimate：a nonparametric retransformation method.J Am Stat Assoc，1983，78（383）：605-610.

［35］ WG Manning JM.Estimating log models：to transform or not to transform.J Health Econ，2001，20（4）：461-494.

［36］Mullahy J.Much ado about two：reconsidering retransformation and the two-part model in health econometrics.J Health Econ，1998，17（3）：247-281.

［37］WG Manning CNM.A two-part model of the demand for medical care：preliminary results from the health insurance study.Health，economics，and health economics，1981：103-123.

［38］L Liu RLS，Johnson BA.Analyzing repeated measures semi-continuous data，with application to an alcohol dependence study.Stat Methods Med Res，2012：0962280212443324.

［39］L Liu JZM，Johnson BA.A multi-level two-part random effects model，with application to an alcohol-dependence study.Stat Med，2008，27（18）：3528-3539.

［40］Su L，Brown S，Ghosh P，et al.Modelling Household Debt and Financial Assets：A Bayesian Approach to a Bivariate Two-Part Model.The University of Sheffield，Department of Economics，2012.

［41］A Zeileis CK，Jackman S.Regression models for count data in R.Journal of Statistical Software，2008，27（8）.［42］BH Neelon AJO，Normand S.A Bayesian model for repeated measures zero-inflated count data with application to outpatient psychiatric service use.Stat Modelling，2010，10（4）：421-439.

［43］P Ghosh PSA.A Bayesian analysis for longitudinal semicontinuous data with an application to an acupuncture clinical trial.Computational statistics＆data analysis，2009，53（3）：699-706.

*：國家自然科學(xué)基金（81172741；30972537）

△通信作者：劉艷，E-mail：liuyan＠ems.hrbmu.edu.cn

鄧 妍）