零膨脹計數(shù)數(shù)據(jù)回歸模型的選擇與比較及R語言的實現(xiàn)*

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

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

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

本研究擬比較不同的算法模型以及傳統(tǒng)的統(tǒng)計模型，在零膨脹數(shù)據(jù)預(yù)測分析上的優(yōu)劣，從而為零膨脹數(shù)據(jù)的分析提供一個新的思路。

方法與實例分析

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

表1　HIV合并血友病數(shù)據(jù)變量基本描述

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

圖1　deaths變量取值分布

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

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

(1)

(2)

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

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

該函數(shù)基本格式如下：

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

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

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

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

表2　零膨脹Poisson回歸模型結(jié)果

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

(3)

(4)

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

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

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

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

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

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

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

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

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

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

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

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

圖2　不同模型擬合結(jié)果比較

小　　結(jié)

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

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

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