零膨脹計數數據回歸模型的選擇與比較及R語言的實現*

劉振球　嚴　瓊　左佳鷺　方綺雯　張鐵軍△

【提　要】　目的　探討和比較不同模型在零膨脹數據回歸分析中的應用。 方法　在R語言中，擬合HIV合并血友病數據的對數線性模型、零膨脹模型、隨機森林、決策樹以及支持向量機模型，通過比較標準化均方誤差和均方根誤差，對模型進行評價與選擇。結果　從標準化均方誤差和均方根誤差來看，隨機森林是對原始數據擬合的最好的模型，隨后是支持向量機和決策樹模型，而經典的計數模型表現則相對較差。結論　在對零膨脹計數資料進行回歸預測時，機器學習方法的效果優于經典的計數模型。

在醫學研究中，經常遇到因變量為計數資料的問題，如一段時間內癲癇患者的抽搐發作次數、某種化療藥物引起患者的嘔吐次數等。這些事件發生次數取值為非負整數，且多服從正偏態分布。研究此類問題時，我們通常假設事件發生次數服從Poisson分布，將該關注事件發生的頻數作為因變量，擬合廣義線性Poisson回歸模型，去探索其他自變量對該事件的影響[1]。當事件發生次數過度離散時，我們可以選擇負二項回歸模型。但是，在某些罕見事件的研究中，常遇到許多觀察個體在研究時間內，并未發生該結局事件，因此數據中會有相當比例的結局變量取值為零，并且零的比例超過Poisson回歸和負二項回歸的預測能力，故稱為零膨脹[2]。Lambert首次建立了零膨脹Poisson回歸模型[3]。該模型的提出，有效地解決了零膨脹數據的分析問題，使分析結果更加準確。近年，關于零膨脹計數資料的模型選擇問題已成為一個研究熱點，各種模型層出不窮，比如Hurdle計數模型[4]，半連續數據兩部模型等[5]。此類模型均是基于經典的統計方法，雖然能夠用明確的數學表達式對原始數據進行展示和解釋，但是由于受限于原始數據的結構以及對原始數據的一些假設，因此在對未知數據進行預測時，往往效果不好。隨著計算機技術的不斷發展，機器學習方法也被越來越多地運用到實際問題的處理過程中。以往的研究證實，在時間序列數據的預測上，算法模型，比如隨機森林、支持向量機、決策樹等，其預測效果明顯優于傳統的統計模型[6-7]。

本研究擬比較不同的算法模型以及傳統的統計模型，在零膨脹數據預測分析上的優劣，從而為零膨脹數據的分析提供一個新的思路。

方法與實例分析

本文所采用的分析數據來自于美國國家癌癥研究所資助的多中心血友病隊列研究(http://www.stat.berkeley.edu/users/statlabs/labs.html)。該項研究從1978年1月1日到1995年12月31日在歐美16個治療中心跟蹤隨訪了超過1600個血友病病人，所得數據共有2144個觀測值及6個變量。表1為變量的基本描述。

表1　HIV合并血友病數據變量基本描述

在上述變量中，deaths是一個零膨脹變量，其取值分布如圖1所示。

圖1　deaths變量取值分布

死亡數等于零的組占比為85.5%，因此，該數據為典型的零膨脹數據。在經典統計學的基礎上，我們一般采用零膨脹計數數據模型對其進行回歸。該模型由兩個部分構成，一部分為集中在零點的點質量，如logistic或者probit回歸模型；第二部分為某種計數分布，比如Poisson分布或者負二項分布。以零膨脹Poisson模型為例，其密度函數可以表示為：

p(yi=0|xi)=pi+(1-pi)exp(-μi)

(1)

(2)

上式中，yi為某事件發生數，xi為協變量向量，μi為第i個個體的期望Poisson計數，pi為二項分布產生的零計數概率。零膨脹負二項分布的概率分布與零膨脹Poisson分布模型類似，二者第一部分相同，而在非零部分選用了負二項分布。

R語言pscl包中的zeroinfl()函數可用來擬合零膨脹負二項分布模型、零膨脹Poisson模型，以及零膨脹幾何分布模型。

該函數基本格式如下：

fit_zero <- zeroinfl(deaths～hiv+factor+py+age | hiv+py+age,data,dist)

管道符“|”將模型分為兩個部分，前面是零部分，后面是非零部分，至于每一部分用什么變量進行擬合，無法先驗確定，可以通過多次嘗試來決定。該函數默認是進行零膨脹Poisson回歸，我們可以根據dist參數選擇相應的非零部分的分布模型。

關于選擇何種分布模型，我們可以使用過度離散檢驗[8]和Vuong檢驗[9]來決定，與之對應的函數分別是odtest()和vuong()。除此之外，我們也可以利用AIC，BIC，標準均方誤差以及均方根誤差等統計指標作為判斷標準。

零膨脹Poisson回歸模型得出的結果見表2。

表2　零膨脹Poisson回歸模型結果

為了比較不同的模型，包括經典的計數模型以及機器學習模型，對于該數據的擬合情況，采用標準均方誤差(NMSE)以及均方根誤差(RMSE)作為判斷標準，對各個模型進行評價。NMSE和RMSE的計算公式如下：

(3)

(4)

在模型擬合的過程中，為了保證結果的穩健性，采用了10重交叉驗證，最后對NMSE和RMSE取均值。實現這一過程的函數見附件。該自定義函數命名為zero_fl，包含9個參數，分別為data，model，formula，id，tar，z=10，p=0.8，dist=′poisson′，seed=2017；其含義分別為：

(1)data：傳入的數據集。

(2)model：選擇的模型，接受一個字符串，比如”rf”，表示進行隨機森林擬合。

(3)formula：針對不同模型的回歸公式。

(4)id：接受一個正整數，表示根據這個變量對原始數據集進行均衡切分。

(5)tar：接受一個正整數，表示目的變量。

(6)z：默認值為10，表示進行10重交叉驗證。

(7)p：默認值是0.8，表示進行交叉驗證時，將80%的數據設置為訓練集。

(8)dist：表示進行零膨脹模型時，采用何種分布，可選“poisson”，“bioneg”和“geometric”。

(9)seed：默認是2017，用于設置隨機數種子。

該函數最終返回的是NMSE和RMSE。對于本文中使用的數據集，最終不同模型擬合的結果如圖2和表3所示。

Liner:Poisson對數線性模型；RF：隨機森林；SVM：支持向量機；Rpart：決策樹模型；Poisson：Poisson零膨脹模型；Negbin：負二項分布零膨脹模型；Geometric：幾何分布零膨脹模型

圖2　不同模型擬合結果比較

小　　結

計數數據是我們在醫學科研中經常遇到的一個問題，對于此類問題，常用的方法是廣義Poisson對數線性模型[10]。但是對于因變量中零過多的情況，傳統的統計模型則不能對數據進行很好的擬合，從而造成數據信息使用不全，導致偏離甚至錯誤的結論。零膨脹模型的提出很好地解決了這個問題，我們可以根據原始數據的特征，比如零在結局變量中所占的比例，選擇相應的零膨脹模型；也可以對多個不同的零膨脹模型進行統計學比較，從而選擇最優模型。

從上文中的結果來看，專門為計數資料設計的若干經典統計模型的表現整體不如廣譜的算法模型。擬合程度最好的是隨機森林，隨后是支持向量機和決策樹模型，表現最差的是Poisson對數線性模型，零膨脹模型介于此二類模型之間。這提示我們今后在處理零膨脹數據時，如果需要對每一個自變量進行解釋，則可以選擇合適的零膨脹回歸模型，如果需要對未知數據進行預測，機器學習方法是一個更好的選擇。

經典的統計模型與算法模型存在本質區別，前者通常要求數據滿足若干假設，如果數據滿足所有假設，則經典模型會擬合出完美的結果，數學上也能被精確描述，我們從而能夠根據模型對數據和結果作出正確的解釋。而機器學習算法不基于原始數據的任何假設，因此適用范圍更加廣泛。這些方法預測精度高，但是不會得到類似P值那樣的顯著性度量指標，也無法用精確的數學公式來描述，更不會用諸如無偏性等概念來評價模型，所以交叉驗證的方法被廣泛用于評價算法模型。本文中提供的R語言代碼，簡單實現了對不同模型的10重交叉驗證，有助于我們快速得到準確的結果。