解決交通事故數(shù)據(jù)分析中零值問題的模型

徐 建，孫 璐

(1.東南大學(xué) 交通學(xué)院，南京210096;2.美國德克薩斯州大學(xué) 奧斯汀分校交通研究中心，奧斯汀78712;3.美國華盛頓Catholic 大學(xué) 土木工程系，華盛頓20064)

0 引 言

應(yīng)用計量模型回歸分析交通事故已成為交通安全領(lǐng)域的重要研究內(nèi)容。由于交通事故具有離散、獨立等特性，泊松模型最早被用于研究，后來衍生出大量回歸模型。人們在研究回歸模型時，主要研究模型中隨機(jī)項的選取問題以提高擬合效果，卻往往忽視對數(shù)據(jù)集本身的分析。由于交通事故(特別是死亡事故)具有零星、隨機(jī)、小概率等特征，研究單元(如人口普查區(qū)、路段等)在一定時期內(nèi)往往存在大量零值現(xiàn)象。常規(guī)的回歸模型如泊松模型、泊松伽馬模型或負(fù)二項模型等都是假設(shè)交通事故發(fā)生次數(shù)服從伯努利分布，但由于外界因素如環(huán)境、地形等，導(dǎo)致某些路段從未發(fā)生事故，這與原假設(shè)產(chǎn)生矛盾，因此用常規(guī)模型來研究大量零值問題，往往無法較好地擬合，甚至?xí)a(chǎn)生錯誤的結(jié)果。

近些年，為了解決這種大量零值問題，研究人員提出了多種回歸模型，主要包括3 類，分別是零膨脹模型、多層零膨脹模型和Lindley 模型，每類模型又涉及泊松和負(fù)二項兩種主要分布。①零膨脹模型。該模型分為兩個過程，一個過程假設(shè)只產(chǎn)生零值，即觀察值的數(shù)量為零，另一個過程假設(shè)觀察值與泊松或負(fù)二項分布一致。該模型由Lambert 于1992 年首次提出［1］，并于近些年得到廣泛應(yīng)用，如卡車事故分析［2］、摩托車事故分析［3］、行人事故分析［4］、兩車道路段上車輛事故分析［5］、死亡事故分析［6］等。而國內(nèi)學(xué)者較少涉及，其中馬壯林等運用零膨脹模型分析了交通事故起數(shù)與時間、道路及交通環(huán)境等潛在影響因素之間的關(guān)系，從時空角度，構(gòu)建交通事故起數(shù)時段、周日和月分布模型［7］。②多層零膨脹模型。多層零膨脹模型主要應(yīng)用于離散、分層且含有大量零值的數(shù)據(jù)集。目前該模型多應(yīng)用于生物醫(yī)學(xué)等領(lǐng)域［8-9］，鮮有用于交通安全領(lǐng)域，特別是交通事故分析。該模型主要包括多層零膨脹泊松和多層零膨脹負(fù)二項兩種模型，模型擬合過程與零膨脹模型類似，也分為兩個過程。③Lindley 模型。該模型由Lindley 于1958 年提出［10］，但直到最近幾年才由Zamani 等學(xué)者［11-13］應(yīng)用于交通事故分析。不同于零膨脹模型和多層零膨脹模型的兩個擬合過程，Lindley 模型只有一個過程，該過程將Lindley 分布融入到泊松模型或負(fù)二項模型。

目前，國內(nèi)外鮮有文獻(xiàn)對這3 類6 種模型的綜合分析，研究適用條件和比較擬合效果。基于此，本文首先系統(tǒng)闡述這3 類6 種模型的結(jié)構(gòu)特征、適用條件、參數(shù)估計等。其次介紹多種擬合優(yōu)度措施，包括離差信息準(zhǔn)則(DIC)、平均絕對偏差(MAD)、預(yù)測均方誤差(MSPE)和最大累計殘差(MCPD)以及交叉驗證評價法(CV)。最后以美國某州2010 年交通事故為實例綜合比較和分析。

1 模型設(shè)定

為了敘述方便，模型名稱以英文字母替代，Poisson 模型為泊松模型;NB 模型為負(fù)二項模型;ZIP 模型為零膨脹泊松模型;ZINB 模型為零膨脹負(fù)二項模型;MZIP 模型為多層零膨脹泊松模型;MZINB 模型為多層零膨脹負(fù)二項模型;Poisson-L模型為泊松Lindley 模型;NB-L 為負(fù)二項Lindley模型。

最早應(yīng)用于交通事故研究的回歸模型是Poisson 模型:

式中:Yi為隨機(jī)變量，表示一定時期內(nèi)路段i 上發(fā)生的交通事故數(shù)，i=1，2，3，…，n，即共有n 條路段;yi為Yi的具體觀察值;均值與方差相等，即E(Yi)=Var(Yi)=λi。

那么，交通事故數(shù)與協(xié)變量(即交通事故影響因素)之間的關(guān)系如式(2)所示:

式中:β0為截距項，即常數(shù)項;βk為第k 個協(xié)變量的系數(shù);xik為對應(yīng)于第i 條路段的第k 個協(xié)變量。

近些年研究人員發(fā)現(xiàn)交通事故與某個(或某些)協(xié)變量存在主要關(guān)系，同時由該變量可以構(gòu)建其他變量，即與其他變量存在一定的相關(guān)性，為此提出設(shè)置曝光變量。本文采用VMT 即vehiclemiles traveled(車輛路段行駛距離)作為曝光變量。另外，考慮到數(shù)據(jù)錯誤、計算誤差等，本文引入誤差項，如式(3)所示:

式中:VMTi為曝光變量，表示第i 條路段上車輛行駛總長度;α 為VMTi的系數(shù);εi為誤差項，服從伽馬分布，即εi～gamma(γ，γ)。

由于Poisson 模型響應(yīng)變量的均值和方差相等，而實際交通事故的方差往往大于均值，所以該模型有較大的局限性。在Poisson 模型的基礎(chǔ)上，NB 模型由于不受這樣的約束，得到廣泛的應(yīng)用，如式(4)所示［14］:

式中:ρ 為過度散布系數(shù)，表示變量的散布程度，當(dāng)ρ=0 時，NB 模型變?yōu)镻oisson 模型。NB 模型變量的均值E(Yi)=λi，方差Var(Yi)=λi+ρλi2。

1.1 零膨脹模型

零膨脹模型［1］的目的是為了解決計數(shù)模型中存在大量零值問題。零膨脹模型分為ZIP 模型和ZINB 模型，模型分為兩個過程，一個過程假設(shè)只產(chǎn)生零值，另外一個過程假設(shè)產(chǎn)生的交通事故數(shù)與Poisson 模型或負(fù)二項模型一致。本文首先介紹ZIP 模型，如式(5)所示:

式中:yi=0 表示ZIP 模型的第一個過程(稱為結(jié)構(gòu)零值)，只產(chǎn)生零值，一般采用對數(shù)分布(logit);yi=1，2，…表示第二個過程(稱為取樣零值)，產(chǎn)生的交通事故數(shù)與Poisson 模型一致。

當(dāng)θi＞0 時，Yi的邊緣分布顯示過度散布。對應(yīng)于第一個過程，即結(jié)構(gòu)零值的均值函數(shù)為對數(shù)線性關(guān)系式，如式(6)所示:

如前面描述，ZIP 模型考慮了過多零值問題，而由于Poisson 分布的局限，未考慮到過度散布問題，而NB 模型則滿足過度散布性質(zhì)，因此將第二個過程由NB 分布代替Poisson 分布，ZIP 模型即變?yōu)閆INB 模型，如式(7)所示:

從式(7)可以看出，Poisson 模型、NB 模型、ZIP 模型和ZINB 模型彼此相關(guān)，例如，當(dāng)ρ =0時，ZINB 模型將變?yōu)閆IP 模型;當(dāng)ρ=0 和θi=0時，ZINB 模型將變?yōu)镻oisson 模型。

以上為ZIP 模型和ZINB 模型的設(shè)定過程。針對零膨脹模型的特性，Vuong 在1989 年提出了Vuong 統(tǒng)計檢驗［15］，用于比較ZIP 模型、ZINB 模型與Poisson 和NB 模型的優(yōu)劣。Vuong 統(tǒng)計檢驗可以在統(tǒng)計軟件中直接實現(xiàn)，如STATA 軟件，或通過編程計算。

1.2 多層零膨脹模型

多層零膨脹模型主要應(yīng)用于離散、分層且含有大量零值的數(shù)據(jù)集，目前該模型多應(yīng)用于生物醫(yī)學(xué)等領(lǐng)域，鮮有用于交通事故分析。該模型主要包括MZIP 模型和MZINB 模型，模型分布過程與零膨脹模型類似，共分為兩個過程，一個過程假設(shè)只產(chǎn)生零值，另一個過程假設(shè)產(chǎn)生的交通事故數(shù)與Poisson 模型或NB 模型一致［9］。

首先介紹MZIP 模型。由于ZIP 模型和ZINB模型未考慮分層以及同一層級上的相互關(guān)系(非獨立性)，因此本文引入兩個隨機(jī)項ηlogiti和ηPi，分別代表logit 過程和Poisson 過程中第i 層級上觀察值之間的依賴程度。設(shè)Yijk表示在第i 層級第j 個路段上第k 個觀察值，其中i=1，2，…，m;j=1，2，…，ni;k=1，2，…，nij。假設(shè)不同層級上的觀察值相互獨立，而同一層級上的觀察值存在相互關(guān)系，因此表達(dá)式為:

式中:wijk和xijk分別為logit 和Poisson 過程協(xié)變量;分別是logit 和Poisson 在第j路段第i 層級上的隨機(jī)誤差項:

那么有:

類似于MZIP 模型，MZINB 模型兩個過程的均值函數(shù)表達(dá)式為:

多層零膨脹模型兩個過程的參數(shù)通過EM 算法估計［9］，即模型的對數(shù)似然函數(shù)由固定數(shù)值EM 算法最大化，進(jìn)而求得相應(yīng)的參數(shù)估計。

1.3 Lindley 模型

不同于零膨脹模型和多層零膨脹模型的兩個擬合過程，Lindley 模型只有一個過程，該過程將Lindley 分布融入到Poisson 模型(Poisson-L 模型)或NB 模型(NB-L 模型)，如式(16)和式(17)所示。這種組合分布具有厚尾性，當(dāng)數(shù)據(jù)集含有大量零值時，能較好地擬合［13］。

式中:f(λ;ρ，φ)表示f 是變量λ 的分布，ρ 和φ為參數(shù);f 服從Lindley 分布。Lindley 分布是由指數(shù)分布和伽馬分布組合而成，如式(18)所示:

均值分別為:

假設(shè)交通事故數(shù)Y 服從Poisson-L 模型或NB-L 模型，那么均值方程為:

可以看出，式(22)與前面描述的式(3)一致。模型的方差為:

關(guān)于Poisson-L 模型和NB-L 模型的參數(shù)估計詳見文獻(xiàn)［11］，通過WinBUGS 軟件中用馬爾可夫鏈蒙特卡洛法進(jìn)行參數(shù)估計。

2 模型比較

本文引入多種擬合優(yōu)度，并采用交叉驗證評價法(CV)對模型的擬合效果進(jìn)行比較。

(1)離差信息準(zhǔn)則(DIC)

DIC 基于貝葉斯理論并考慮了模型的復(fù)雜性，它與赤池信息準(zhǔn)則(AIC)意義相同［15］，表達(dá)式為:

(2)平均絕對偏差(MAD)

MAD 是一種模型平均錯誤預(yù)測的測量方法，與平均預(yù)測誤差(MPB)不同，MAD 的值越接近0時，說明模型預(yù)測效果越好［16］，表達(dá)式為:

式中:n 為樣本量大小。

(3)預(yù)測均方誤差(MSPE)

MSPE 與均方誤差類似，是一種衡量觀察事故率與預(yù)測事故率差異的方法，同時可以評價數(shù)據(jù)的變化程度，MSPE 的值越小，說明預(yù)測模型具有更好的精確度，表達(dá)式為:

式中:n2為樣本量大小。

(4)最大累計殘差(MCPD)

MCPD 是偏離0 的累計殘差(觀察事故率與預(yù)測事故率的差值)最大絕對值［17］。MCPD 的值越小，說明模型擬合效果越好。

(5)交叉驗證評價法(CV)

CV 是一種基于貝葉斯理論，用于評價模型復(fù)雜性和擬合情況的方法，它能靈活可靠地檢驗不同模型的適合度，特別針對大量零值問題。為了減少馬爾可夫鏈蒙特卡洛算法的運算時間，本文引入K-折交叉驗證。K-折交叉驗證的預(yù)測能力表達(dá)式可在文獻(xiàn)［18］中查詢，本文不再贅述。

3 實例分析

3.1 數(shù)據(jù)概述

本文以美國某州交通事故為實例，選取了該州2010 年發(fā)生在州級道路系統(tǒng)內(nèi)(由州交通廳管理的路網(wǎng))事故數(shù)據(jù)集，共計243 388 起。按照同性質(zhì)要求(同一路段的限速、交通量、路面寬度、車道數(shù)、平縱線形等一致)，該州117472.28 km 道路網(wǎng)被劃分成277 510 條同性質(zhì)路段，平均路段長度為0.42 km，其中90%以上的路段長度集中于0 ～1.6 km。

需要說明的是:零膨脹模型和Lindley 模型的研究單元均為同性質(zhì)路段，而多層零膨脹模型的研究單元分為兩個層級:第一層級為人口普查區(qū)，第二層級為同性質(zhì)路段。

利用ArcGIS 軟件將243 388 起交通事故與同性質(zhì)路段合并，即可得到每條路段的交通事故數(shù)。將事故數(shù)按數(shù)量進(jìn)行劃分，如表1 所示。從表中可以看到，82.05%的同性質(zhì)路段上事故數(shù)為零，89.36%的同性質(zhì)路段未發(fā)生受傷事故和死亡事故，因而，該數(shù)據(jù)集存在大量零值現(xiàn)象。

表1 路段上交通事故數(shù)分布Table 1 Distributions of crash counts on segments

3.2 因素選取

本文中響應(yīng)變量為交通事故數(shù)，曝光變量為車輛路段行駛距離VMT(VMT=AADT×路段長度×365)，而其他協(xié)變量主要包括幾何線形、交通特性和天氣狀況3 類共8 種，這些協(xié)變量的均值、標(biāo)準(zhǔn)差、最小值和最大值如表2 所示。這里需要說明的是，平均路肩寬度指外側(cè)路肩和內(nèi)側(cè)路肩的平均值;中央分隔帶寬度不包含內(nèi)側(cè)路肩寬度;平曲線度數(shù)指每一百英尺對應(yīng)的曲線度數(shù)。

表2 路段的變量統(tǒng)計結(jié)果Table 2 Summary statistics of variables for segments

3.3 結(jié)果分析

本文運用基于貝葉斯理論的馬爾可夫鏈蒙特卡洛算法估計模型參數(shù)［19］，推演過程共有6 條馬爾可夫鏈，每條鏈迭代20 000 次，前5000 次迭代作為burn-in 過程，以消除初始參數(shù)的影響，余下15 000 次迭代用于參數(shù)估計;同時，本文采用Gelman-Rubin(G-R)收斂統(tǒng)計來驗證收斂適用性，G-R 收斂統(tǒng)計值小于1.1。整個模型估計在WinBUGS 軟件中完成，最終回歸結(jié)果如表3 所示。注:ZIP 模型、ZINB 模型、MZIP 模型和MZINB 模型的logit 過程，即結(jié)構(gòu)零值過程的系數(shù)和標(biāo)準(zhǔn)差未顯示在表格中。

從表3 可知:6 種模型的協(xié)變量系數(shù)與標(biāo)準(zhǔn)差各不相同，特別是有些協(xié)變量的邊際效應(yīng)在不同模型中甚至出現(xiàn)完全相反的結(jié)果，如車道數(shù)和降雨量對交通事故的影響，在ZIP 模型和MZIP 模型中成負(fù)相關(guān)，而在其他4 種模型中則成正相關(guān);另外，車道年平均日交通量在ZIP 模型、ZINB 模型、MZIP 模型和Poisson-L 模型中與交通事故成負(fù)相關(guān)，而在MZINB 和NB-L 模型中成正相關(guān)。這充分說明不同模型會有不同的參數(shù)估計，產(chǎn)生不同的擬合效果。

從表3 同樣可以得到:曝光變量(VMT)和車道數(shù)、車道年平均日交通量以及降雨量等協(xié)變量與交通事故數(shù)成正相關(guān)，即交通事故數(shù)隨著協(xié)變量數(shù)值增大而增加，而其他協(xié)變量成負(fù)相關(guān)。值得注意的是，從直覺來講最大限速越高，事故越多，但統(tǒng)計結(jié)果卻完全相反，即并非限速越高事故越多，可能由于限速高的路段線形往往更好，同時駕駛員注意力更集中，事故反而少。表4 為各模型擬合優(yōu)度比較結(jié)果。

表3 各交通事故模型的估計結(jié)果Table 3 Estimation results of models for crash count

從表4 中可以看出:NB-L 模型的離差信息準(zhǔn)則(DIC)最小，其他依次是MZINB 模型、ZINB 模型、Poisson-L 模型、MZIP 模型，最大是ZIP 模型。由于DIC 越小說明模型擬合效果越好，由此得到，NB-L 模型的擬合效果最好，其次是MZINB 模型。類似于離差信息準(zhǔn)則(DIC)，通過比較平均絕對偏差(MAD)、預(yù)測誤差均方(MSPE)和最大累計殘差(MCPD)，均可得到相同的結(jié)論。需要說明的是，本文采用的數(shù)據(jù)集具有較強(qiáng)的過度散布性(從表2 可以看出，交通事故數(shù)的方差遠(yuǎn)大于均值)，因而NB 模型比Poisson 模型更合適。若將NB 模型和Poisson 模型分開比較，可知Lindley 模型比多層零膨脹模型和零膨脹模型更合適，而多層零膨脹模型又比零膨脹模型更合適。

表4 各交通事故模型擬合優(yōu)度比較Table 4 Comparison of goodness-of-fit of models for crashes

除了比較4 種擬合優(yōu)度(DIC、MAD、MSPE 和MCPD)外，本文還選取了交叉驗證評價法(CV)，圖1 表示各模型的CV 預(yù)測能力，注:由于路段最大事故數(shù)達(dá)到387 起(見表2)，圖1 僅顯示同性質(zhì)路段上事故數(shù)為0 ～10 起，超過10 個事故數(shù)的未顯示。從圖中可以看出，雖然針對不同事故數(shù)，各模型的預(yù)測能力不盡相同，但總體上NB-L 模型的預(yù)測能力最強(qiáng)，MZINB 模型次之，而ZIP 模型預(yù)測能力最弱，這與前面4 種擬合優(yōu)度的比較結(jié)果類似。

4 結(jié)束語

針對交通事故數(shù)據(jù)分析中大量零值問題，本文對目前已有的3 類6 種模型綜合分析，發(fā)現(xiàn)ZIP 模型、ZINB 模型、MZIP 模型和MZINB 模型分為兩個過程(結(jié)構(gòu)零值和取樣零值)，其中MZIP 模型和MZINB 模型適合于分層數(shù)據(jù)集，而Poisson-L 模型和NB-L 模型只有一個過程，這更符合統(tǒng)計學(xué)理論;同時，通過離差信息準(zhǔn)則(DIC)、平均絕對偏差(MAD)、預(yù)測誤差均方(MSPE)和最大累計殘差(MCPD)等4 種擬合優(yōu)度以及交叉驗證評價法(CV)，以美國某州2010 年交通事故為實例，運用馬爾可夫鏈蒙特卡洛算法，結(jié)果表明，3 類模型中Lindley 模型擬合效果最好，多層零膨脹模型其次，而6 種模型中負(fù)二項Lindley 模型擬合效果最好，MZINB 模型其次。本文可為運用計數(shù)模型研究大量零值現(xiàn)象提供借鑒和參考。

圖1 各模型的CV 預(yù)測能力比較Fig.1 Model comparison of predictive abilities using cross-validation

［1］Lambert D.Zero-inflated poisson regression，with an application to defects in manufacturing［J］.Technometrics，1992，34(1):1-14.

［2］Miaou，S.The relationship between truck accidents and geometric design of road sections:poisson versus negative binomial regressions［J］.Accident Analysis＆Prevention，1994，26(4):471-482.

［3］Lee J，Mannering F L.Impact of roadside features on the frequency and severity of run-off-road accidents:an empirical analysis［J］.Accident Analysis＆Prevention，2002，34(2):349-361.

［4］Shankar V N，Ulfarsson G F，Pendyala R M，et al.Modeling crashes involving pedestrians and motorized traffic［J］.Safety Science，2003，41(7):627-640.

［5］Qin X，Ivan J N，Ravishankar N.Selecting exposure measures in crash rate prediction for two-lane highway segments［J］.Accident Analysis＆Prevention，2004，36(2):183-191.

［6］Yan Xue-dong，Wang Bin，An Mei-wu，et al.Distinguishing between rural and urban road segment traffic safety based on zero-inflated negative binomial regression models［DB/OL］.［2013-01-18］.http://www.hindawi.com/journals/ddns/2012/789140/.

［7］馬壯林，邵春福，胡大偉，等.高速公路交通事故起數(shù)時空分析模型［J］.交通運輸工程學(xué)報，2012，12(2):93-99.Ma Zhuang-lin，Shao Chun-fu，Hu Da-wei，et al.Temporal-spatial analysis model of traffic accident frequency on expressway［J］.Journal of Traffic and Transportation Engineering，2012，12(2):93-99.

［8］Lee A H，Wang Kui，Scott J A，et al.Multi-level zeroinflated Poisson regression modeling of correlated count data with excess zeros［J］.Statistical Methods in Medical Research，2006，15(1):47-61.

［9］Moghimbeigi A，Eshraghian M R，Mohammad K，et al.Multilevel zero-inflated negative binomial regression modeling for over-dispersed count data with extra zeros［J］.Journal of Applied Statistics，2008，35(10):1193-1202.

［10］Lindley D V.Fiducial distributions and Bayes＇theorem［J］.Journal of the Royal Statistical Society，1958，20(1):102-107.

［11］Zamani H，Ismail N.Negative binomial-Lindley distribution and its application［J］.Journal of Mathematics and Statistics，2010，6(1):4-9.

［12］Lord D，Geedipally S R.The negative binomial-Lindley distribution as a tool for analyzing crash data characterized by a large amount of zeros［J］.Accident Analysis and Prevention，2011，43(5):1738-1742.

［13］Reddy G S，Lord D，Dhavala S S.The negative binomial-Lindley generalized linear model:Characteristics and application using crash data［J］.Accident Analysis＆Prevention，2012，45(3):258-265.

［14］Miaou S.The relationship between truck accidents and geometric design of road sections:Poisson versus negative binomial regressions［J］.Accident Analysis＆Prevention，1994，26(4):471-482.

［15］Vuong Q H.Likelihood ratio tests for model selection and non-nested hypothesis［J］.Econometrica，1989，57(2)，307-333.

［16］Oh J，Lyon C，Washington S P，et al.Validation of the FHWA crash models for rural intersections:lessons learned［C］∥Transportation Research Record，Transportation Research Board，National Research Council，Washington DC，2003:41-49.

［17］Geedipally S R，Patil S，Lord D.Examination of methods for estimating crash counts according to their collision type［J］.Transportation Research Record，Transportation Research Board，National Research Council，Washington DC，2010:12-20.

［18］Huang H L，Chin H C.Modeling road traffic crashes with zero-inflation and site-specific random effects［J］.Statistical Methods＆Applications，2010，19(3):445-462.

［19］Thomas A，O＇Hara B，Ligges U，et al.Making BUGS open［J］.R News，2006，6(1):12-17.