三種模型對廣東省傷寒副傷寒逐月發(fā)病數(shù)預(yù)測的比較*

李 驪錢 俊楊 軍歐春泉△

三種模型對廣東省傷寒副傷寒逐月發(fā)病數(shù)預(yù)測的比較*

李 驪1錢 俊2楊 軍1歐春泉1△

目的應(yīng)用三種統(tǒng)計(jì)模型預(yù)測傷寒副傷寒的發(fā)病趨勢，比較其預(yù)測效果，為傷寒副傷寒的預(yù)測和防控提供科學(xué)依據(jù)。方法利用廣東省2008年5月至2012年4月四年的傷寒副傷寒逐月發(fā)病資料，分別擬合季節(jié)性綜合自回歸滑動平均（SARIMA）模型、經(jīng)傅里葉季節(jié)性調(diào)整的綜合自回歸滑動平均（FARIMA）模型和動態(tài)諧波回歸（DHR）模型，并用前面建立的三種模型預(yù)測后續(xù)半年（2012年5月—2012年10月）的逐月發(fā)病數(shù)。結(jié)果傷寒副傷寒的發(fā)病有明顯的周期性和季節(jié)特征，周期為1年，7-8月份為發(fā)病高峰期。流行強(qiáng)度和流行高峰出現(xiàn)的月份均存在一定的年度差異。三種模型擬合四年的傷寒副傷寒發(fā)病情況，其平均絕對百分比誤差（MAPE）依次為：DHR模型（7.8%）＜FARIMA模型（12.9%）＜SARIMA模型（13.4%）；三種模型預(yù)測后續(xù)半年的發(fā)病情況，其MAPE依次為：DHR模型（3.5%）＜FARIMA模型（5.6%）＜SARIMA模型（6.8%），其他模型評價(jià)指標(biāo)結(jié)果也類似。結(jié)論三種方法均有較佳的預(yù)測效果。相對而言，DHR的預(yù)測精度更高。本研究可為常見傳染病的預(yù)測提供一定的方法學(xué)參考。

SARIMA模型 FARIMA模型 DHR模型 預(yù)測 傷寒副傷寒

傷寒副傷寒是我國法定的乙類傳染病。在世界范圍內(nèi)，該病是導(dǎo)致死亡的重要原因［1］。本病亦是我國常見的腸道傳染病之一，在我國終年散發(fā)，部分地區(qū)出現(xiàn)暴發(fā)流行［2-3］。其傳染性強(qiáng)、病程長、易復(fù)發(fā)、并發(fā)癥多、疾病負(fù)擔(dān)較重［4］。由于流動人口增加、耐藥菌株增多等問題的出現(xiàn)使其防控難度進(jìn)一步增加，因而有效地預(yù)測傷寒副傷寒的發(fā)病情況，對于衛(wèi)生管理部門及時(shí)做好預(yù)防干預(yù)措施具有重要意義。SARIMA模型是分析時(shí)間序列的經(jīng)典時(shí)間域方法。FARIMA模型和DHR模型則是時(shí)間域與頻率域相結(jié)合的方法，二者能結(jié)合時(shí)間域與頻率域兩類時(shí)間序列分析方法的優(yōu)勢，國內(nèi)尚未有研究使用此兩種方法預(yù)測傳染病的發(fā)病。筆者試構(gòu)建以廣東省2008年5月至2012年4月傷寒副傷寒逐月發(fā)病數(shù)據(jù)為基礎(chǔ)的FARIMA模型和DHR模型，并將二者與常用的SARIMA模型進(jìn)行比較，探討傷寒副傷寒預(yù)測的優(yōu)化模型，為建立該病的預(yù)測系統(tǒng)提供科學(xué)依據(jù)，同時(shí)也為其他傳染病的預(yù)測提供思路。

資料與方法

1.資料來源

從廣東省衛(wèi)生廳官方網(wǎng)站獲得全省2008年5月至2012年10月期間傷寒副傷寒逐月發(fā)病資料。廣東省有很好的傳染病監(jiān)測系統(tǒng)，保證了本研究資料的可靠性。

2.統(tǒng)計(jì)方法及建模過程

基于2008年5月至2012年4月傷寒副傷寒逐月發(fā)病數(shù)據(jù)建立以下三種模型，用2012年5月至2012年10月的數(shù)據(jù)對模型的預(yù)測效果進(jìn)行評價(jià)。

（1）SARIMA模型［5］

SARIMA模型的構(gòu)建分以下幾個步驟實(shí)現(xiàn)：

①對數(shù)據(jù)進(jìn)行預(yù)處理，即：xt＝lnyt，其中yt為待分析的時(shí)間序列。

②模型識別，即根據(jù)自相關(guān)（ACF）圖，偏自相關(guān)（PACF）圖及AIC最小化準(zhǔn)則確定模型。

③參數(shù)估計(jì)，并檢驗(yàn)參數(shù)是否有統(tǒng)計(jì)學(xué)意義。

④模型診斷及預(yù)測，通過判斷殘差序列是否為白噪聲和存在自相關(guān)性評價(jià)模型是否合適。用確定的模型進(jìn)行預(yù)測。

（2）FARIMA模型

傅里葉方法的主要原理是用一系列正弦或余弦波近似擬合原始序列的波譜，本例以此擬合序列的季節(jié)組分，即：

其中period表示序列周期，本例由于傷寒副傷寒發(fā)病數(shù)呈現(xiàn)明顯的12個月周期，故設(shè)為12；t表示第t月，取值為1，2，…，n，n表示待分析時(shí)間序列的長度；terms通常取小于等于period/2的整數(shù)。對數(shù)據(jù)進(jìn)行季節(jié)性調(diào)整后構(gòu)建ARIMA模型。

（3）DHR模型［6-7］

動態(tài)諧波回歸（dynamic harmonic regression，DHR）模型是單變量UC（unobserved component）模型的特例，本研究采用的DHR模型為：

其中xt是經(jīng)處理的時(shí)間序列，Tt代表趨勢組分或低頻組分，St表示季節(jié)性周期組分。et為殘差，經(jīng)常將其假設(shè)為服從均數(shù)為0，方差為σ2的正態(tài)分布。由于Tt和St不能直接由原始數(shù)據(jù)得到，故將該模型稱為UC模型。UC模型的各組分可由狀態(tài)空間模型來表達(dá)。在DHR建模中，我們將用到廣義隨機(jī)游走（GRW）類模型。最一般的GRW類模型定義為：

其中，x1，t，x2，t分別是各組分的變化水平和斜率。η1，t，η2，t是對角協(xié)方差陣為常量的正態(tài)隨機(jī)變量。GRW的主要取值方式有隨機(jī)游走（RW）模型：α＝1，β＝γ＝0，η2，t＝0；整合隨機(jī)游走（IRW）模型：0＜α＜1，β＝γ＝0，η2，t＝0等。

在DHR模型中，各組分可由時(shí)變參數(shù)（TVP′s）表示，這些參數(shù)可以采用最優(yōu)卡爾曼濾波（KF）和固定間隔平滑（FIS）方法估計(jì)。在估計(jì)TVP′s時(shí)，需先確定諸如噪聲方差比例（NVR）這些超參數(shù)，超參數(shù)則采用頻率域的方法進(jìn)行估計(jì)。

本研究DHR建模分以下幾個步驟實(shí)現(xiàn)：

①對數(shù)據(jù)進(jìn)行預(yù)處理，即：xt＝lnyt，其中yt為待分析的時(shí)間序列。

②估計(jì)處理后序列的經(jīng)驗(yàn)譜。根據(jù)AIC最小化準(zhǔn)則，自回歸譜階數(shù)設(shè)為24。

③通過分析自回歸譜，確定諧波頻率。

④估計(jì)噪聲方差比例及超參數(shù) IRW模型適用于表達(dá)緩慢變化的情況，RW模型適用于表達(dá)變化快速的情況［6］。故Tt和St的擬合分別采用IRW和RW模型。

⑤利用第④步估計(jì)得到NVR和其他超參數(shù)，由最優(yōu)卡爾曼濾波和固定間隔平滑法實(shí)現(xiàn)每個組分的估計(jì)及預(yù)測功能。

（4）模型評價(jià)指標(biāo)

采用以下指標(biāo)評價(jià)模型的預(yù)測準(zhǔn)確性：均方根誤差（RMSE），平均絕對誤差（MAE），平均絕對百分比誤差（MAPE）及平均絕對偏差百分比（PMAD）。各指標(biāo)的計(jì)算方法如下：

其中，N表示待分析時(shí)間序列的長度；Et表示擬合或預(yù)測誤差，即實(shí)測值與預(yù)測值之差；Yt為時(shí)間序列的實(shí)測值。

（5）統(tǒng)計(jì)軟件

SARIMA和FARIMA建模采用R 2.15.1，DHR建模采用MATLAB 7.6編程實(shí)現(xiàn)。

結(jié) 果

1.發(fā)病趨勢分析

2008-2012年廣東省傷寒副傷寒每月發(fā)病數(shù)見圖1。由圖可以看出，傷寒副傷寒的發(fā)病有明顯的周期性和季節(jié)特征，周期為1年，7～8月份為發(fā)病高峰期，而其他月份也存在一定的病例數(shù)。流行強(qiáng)度和流行高峰出現(xiàn)的月份均存在一定的年度差異，2008-2009年高峰出現(xiàn)在7月，2010-2012年流行高峰在8月份。2010年的發(fā)病人數(shù)最多，季節(jié)性高峰最明顯，全年累計(jì)發(fā)病數(shù)為1828例，2011年又下降至2008年的水平，2012年前三季度總的發(fā)病數(shù)與2011年基本持平，但1～3月的發(fā)病數(shù)較2011年高。

圖1 傷寒副傷寒發(fā)病數(shù)時(shí)序圖

2.各模型構(gòu)建結(jié)果

（1）SARIMA模型

對處理過的序列進(jìn)行單位根檢驗(yàn)，可知序列不平穩(wěn)。進(jìn)一步對序列xt進(jìn)行一次季節(jié)性差分后序列平穩(wěn)。觀察SARIMA（0，0，0）×（0，1，0）12殘差序列的ACF和PACF圖（圖2），提示宜采用乘積混合效應(yīng)模型SARIMA（1，0，0）×（0，1，0）12。模型殘差序列的ACF和PACF圖（圖2）提示模型已經(jīng)很好地消除了序列的自相關(guān)特性，能較好地?cái)M合該時(shí)間序列。模型參數(shù)為：φ1＝-0．504，擬合模型如下：

對系數(shù)進(jìn)行t檢驗(yàn)，得t＝3.603（P＝0.000），說明該系數(shù)有統(tǒng)計(jì)學(xué)意義。模型殘差通過了Ljung-Box檢驗(yàn)，說明殘差為白噪聲。將上式展開，得：

（2）FARIMA模型

首先建立只包含傅里葉項(xiàng)的模型，本例terms取1。

圖2 SARIMA（0，0，0）×（0，1，0）12殘差和SARIMA（1，0，0）×（0，1，0）12殘差的自相關(guān)和偏自相關(guān)圖

對模型的殘差進(jìn)行單位根檢驗(yàn)，可知?dú)埐钚蛄胁黄椒€(wěn)，對殘差進(jìn)行一階差分后序列平穩(wěn)。觀察FARIMA（0，1，0）殘差序列的ACF和PACF圖并依據(jù)AIC最小化準(zhǔn)則選定最終模型FARIMA（0，1，1）。模型參數(shù)為：θ1＝-0．568，a1＝46．419，b1＝-15．676，即得到模型：

對模型系數(shù)進(jìn)行t檢驗(yàn)，三者t值分別為：t＝-3.611（P＝0.000），t＝10.053（P＝0.000），t＝-3.587（P＝0.000），系數(shù)有統(tǒng)計(jì)學(xué)意義。模型殘差通過了Ljung-Box檢驗(yàn)，說明殘差為白噪聲。將上式展開，得到最終模型

（3）DHR模型

模型殘差序列的ACF和PACF圖提示模型已經(jīng)消除了序列的自相關(guān)特性，殘差序列通過了Ljung-Box檢驗(yàn)，說明殘差為白噪聲。綜上，模型能較好地?cái)M合該時(shí)間序列。

3.擬合預(yù)測效果評價(jià)

從四年歷史數(shù)據(jù)的擬合情況來看，三種模型均有較好的擬合效果，DHR模型尤為突出，平均絕對百分比誤差最小（表1）。實(shí)際工作中，準(zhǔn)確預(yù)測傳染病流行高峰發(fā)生的時(shí)間和強(qiáng)度具有更大的公共衛(wèi)生意義。由圖3可見：實(shí)際數(shù)據(jù)在2009年出現(xiàn)雙峰，其他年份均為單峰。SARIMA模型由于受前期數(shù)據(jù)的變化特征影響較大，其擬合值繼2009年后在2010年仍然出現(xiàn)了雙峰。FARIMA的擬合曲線比較平滑，對波峰和波谷的擬合欠佳，特別是2009年擬合的波峰滯后實(shí)際流行高峰一個月，而DHR模型對流行高峰出現(xiàn)的時(shí)間和強(qiáng)度均有較佳擬合效果（圖3）。

用半年的數(shù)據(jù)進(jìn)行前瞻性考核，同樣反映DHR模型有最好的預(yù)測能力，2012年5月至10月的預(yù)測相對誤差均控制在10%以內(nèi)（表2），平均絕對百分比誤差僅有3.5%。FARIMA和SARIMA模型的平均絕對百分比誤差略大，分別為5.6%和6.8%。其他評價(jià)模型的指標(biāo)也得出一致的結(jié)論（表1）。

圖3 三種方法的預(yù)測結(jié)果圖

表1 三種模型的擬合及預(yù)測效果比較

表2 三種模型的預(yù)測結(jié)果

討 論

廣東省傷寒副傷寒高發(fā)于夏秋季，發(fā)病數(shù)在7-8月達(dá)到最高峰，這與江西、云南、貴州、天津等地報(bào)道的情況類似［8-11］，而印度的報(bào)道指出高峰出現(xiàn)在4～6月［12］。本病夏秋季高發(fā)可能是因?yàn)樘鞖鉄?，人們喜歡吃生冷的食物，而這些食物容易受到污染造成人的感染。另外，廣東省夏秋季降水較多，強(qiáng)降雨等極端氣候事件也會增加人群傷寒副傷寒的發(fā)病風(fēng)險(xiǎn)。2012年1～3月的發(fā)病數(shù)較2011年多，可能與當(dāng)年廣東省出現(xiàn)“回南天”有關(guān)。

在短短幾年的研究時(shí)間里，一個省的總?cè)丝跀?shù)通常變化較小，故往往采用逐月發(fā)病人數(shù)的時(shí)間序列數(shù)據(jù)來反映傳染病發(fā)病率的變化趨勢。近年來，國內(nèi)有學(xué)者應(yīng)用指數(shù)曲線和灰色模型預(yù)測傷寒副傷寒疫情的年度變化趨勢［13-14］，因而沒有考慮到時(shí)序的季節(jié)性，而傷寒副傷寒發(fā)病率存在明顯的季節(jié)特征，在進(jìn)行逐月發(fā)病情況的預(yù)測時(shí)不容忽視。另外，有人提出用加權(quán)馬爾可夫鏈預(yù)測傷寒副傷寒的發(fā)病情況［15］，但該法只判斷發(fā)病數(shù)的大致范圍，沒有給出具體的估計(jì)值。本研究分別采用SARIMA模型，F(xiàn)ARIMA模型及DHR模型預(yù)測廣東省傷寒副傷寒的發(fā)病情況。對三種模型的比較發(fā)現(xiàn)：SARIMA模型的擬合效果與預(yù)測精度較低，受前一個周期的季節(jié)性變化規(guī)律影響較大，例如：2009年傷寒副傷寒出現(xiàn)雙峰導(dǎo)致2010年也錯誤地模擬出雙峰。故當(dāng)每個周期的季節(jié)性差異較大時(shí)，不建議使用SARIMA模型。在構(gòu)建FARIMA模型時(shí)涉及諧波的選擇，當(dāng)待選諧波的個數(shù)比較多時(shí)，模型遴選較復(fù)雜，耗時(shí)較長。FARIMA模型假定序列的季節(jié)性特征不隨時(shí)間變化而變化，故當(dāng)各年度季節(jié)特征變化較大時(shí)，擬合效果略遜。DHR模型是一種動態(tài)擬合趨勢組分和季節(jié)性周期組分的方法，它允許時(shí)間序列的季節(jié)性特征隨時(shí)間變化而變化［6］，故能很好地反映序列的變化特征。在本研究中，該法的擬合和預(yù)測效果最好。DHR模型的缺點(diǎn)在于，諧波的選擇比較主觀。有時(shí)根據(jù)AIC準(zhǔn)則確定的自回歸譜并不能準(zhǔn)確地提供有意義的諧波頻率信息，故在使用該法時(shí)，不能盲目地相信單個準(zhǔn)則提供的信息，而應(yīng)該結(jié)合相關(guān)專業(yè)知識和擬合結(jié)果，選擇合理的參數(shù)。一般來說，傳染病發(fā)病資料會呈現(xiàn)一定的季節(jié)性，當(dāng)數(shù)據(jù)既包含長期趨勢又存在季節(jié)性時(shí)，可考慮采用DHR模型對傳染病的發(fā)病情況進(jìn)行短期預(yù)測。當(dāng)季節(jié)性不變或變化不大時(shí)，也可嘗試用FARIMA模型和SARIMA模型建模。本文的預(yù)測都是基于發(fā)病數(shù)本身的，并沒有考慮其他因素對發(fā)病情況的影響。事實(shí)上，傳染病的發(fā)生或多或少與氣象、環(huán)境等因素有一定的聯(lián)系，有待將來進(jìn)一步研究。也可以考慮建立一套指標(biāo)體系，并通過綜合運(yùn)用指標(biāo)體系的方法對傳染病的情況進(jìn)行分析和評價(jià)，確認(rèn)發(fā)生危險(xiǎn)的可能性及嚴(yán)重程度，決定是否發(fā)出危機(jī)報(bào)警［16］。

我國已經(jīng)形成了完善的傳染病監(jiān)測報(bào)告體系，各級衛(wèi)生管理部門及時(shí)發(fā)布統(tǒng)計(jì)數(shù)據(jù)，為預(yù)測模型的建立與修正提供了數(shù)據(jù)支持。本研究采用三種模型預(yù)測廣東省傷寒副傷寒的發(fā)病數(shù)，分析比較了三種模型的優(yōu)劣，指出DHR模型具有更好的預(yù)測準(zhǔn)確性，為傷寒副傷寒發(fā)病情況的預(yù)報(bào)提供了科學(xué)依據(jù)，同時(shí)也為其他常見傳染病的預(yù)測提供方法學(xué)的思路和借鑒。

1.Crump JA，M intz ED.Global trends in typhoid and paratyphoid Fever. Clin Infect Dis，2010，50（2）：241-246.

2.天津醫(yī)學(xué)院流行病學(xué)教研室主編.常見腸道傳染病的防治.天津：天津人民出版社，1973：1-2.

3.衛(wèi)生部疾病預(yù)防控制局，中國疾病預(yù)防控制中心主編.傷寒、副傷寒防治手冊.北京：人民衛(wèi)生出版社，2006：1-2.

4.王魯茜，闞飆.傷寒、副傷寒的全球流行概況及其預(yù)防控制.疾病監(jiān)測，2007，22（7）：492-494.

5.彭志行，陶紅，賈成梅，等.時(shí)間序列分析在麻疹疫情預(yù)測預(yù)警中的應(yīng)用研究.中國衛(wèi)生統(tǒng)計(jì)，2010，27（5）：459-463.

6.Jiang B，Liang SL，Wang JD，et al.Modeling MODIS LAI time series using three statistical methods.Remote Sens Environ，2010，114（7）：1432-1444.

7.Young PC，Pedregal DJ，Tych W.Dynam ic harmonic Regression.Journal of Forecasting，1999，18（6）：369-394.

8.劉曉青，袁輝，王健，等.1998-2007年江西省傷寒副傷寒流行病學(xué)分析.海峽預(yù)防醫(yī)學(xué)雜志，2009，15（2）：27-28.

9.黃玉芬，劉志濤.云南省2008-2010年傷寒副傷寒發(fā)病情況.職業(yè)與健康，2011，27（18）：2131-2132.

10.黎明，唐光鵬，姚光海，等.貴州省2004-2007年腸道傳染病流行特征分析.醫(yī)學(xué)動物防制，2012，28（1）：5-8.

11.吳偉慎，解曉華，單愛蘭，等.天津市1990-2005年傷寒副傷寒流行趨勢分析.現(xiàn)代預(yù)防醫(yī)學(xué)，2007，34（17）：3304-3305.

12.Kanungo S，Dutta S，Sur D.Epidem iology of typhoid and paratyphoid fever in India.J Infect Dev Ctries.2008，2（6）：454-60.

13.春雅麗，黃中敏，錢文平.應(yīng)用指數(shù)曲線預(yù)測傷寒副傷寒疫情.上海預(yù)防醫(yī)學(xué)雜志，2005，17（10）：466-467.

14.黎景雪，王培承，邱瑞香，等.灰色模型在我國傷寒副傷寒發(fā)病率預(yù)測中的應(yīng)用.數(shù)理醫(yī)藥學(xué)雜志，2010，23（5）：506-508.

15.彭志行，鮑昌俊，趙楊，等.加權(quán)馬爾可夫鏈在傷寒副傷寒發(fā)病情況預(yù)測分析中的應(yīng)用.中國衛(wèi)生統(tǒng)計(jì)，2008，25（3）：226-229.

16.胡世雄，邢慧嫻，鄧志紅.我國傳染病的預(yù)測預(yù)警現(xiàn)狀.中華預(yù)防醫(yī)學(xué)雜志，2007，41（5）：407-410.

（責(zé)任編輯：劉 壯）

Com parison of Three M odels on Prediction of Incidence of Typhoid Fever and Paratyphoid Fever

Li Li，Qian Jun，Yang Jun，et al（Department of Biostatistics，School of Public Health and Tropical Medicine，Southern Medical University（510515），Guangzhou）

ObjectiveTo compare the performance of three statisticalmethods in forecasting the incidence of typhoid fever and paratyphoid fever.MethodsUsing monthly data of cases w ith typhoid fever and paratyphoid fever in Guangdong Province from May 2008 to April2012，we fitted threemodels，including Seasonal Autoregressive Integrated Moving Average（SARIMA）model，ARIMA model applied to the seasonally adjusted data from Fourier terms（FARIMA）and Dynamic Harmonic Regression（DHR）model.Afterwards，we used the data from May 2012 to October 2012 to assess the accuracy of prediction for these threemodels.ResultsThe incidence of typhoid fever and paratyphoid fever is characterized by an apparent cyclic pattern w ith a one-year seasonal cycle.Themaximum number of cases occurred during July to August.The epidem ic strength and peak differed by years.Mean Absolute Percentage Error（MAPE）of the four-year time series fitted w ith threemethodswas：DHR（7.8%）＜FARIMA（12.9%）＜SARIMA（13.4%），and MAPE of the halfa-year time series forecasted had the same trend：DHR（3.5%）＜FARIMA（5.6%）＜SARIMA（6.8%）.Other index for forecasting accuracy conveyed the sim ilarmessage.ConclusionA ll threemodels had satisfactory prediction capacity.DHRmodel is superior to FARIMA model and SARIMA modelw ith a higher forecast precision.Themethods and findings in this study may throw light on predicting the incidence of other infectious diseases.

SARIMA model；FARIMA model；DHRmodel；Forecasting；Typhoid fever and Paratyphoid fever

*：國家自然科學(xué)基金（81102207）；廣東省自然科學(xué)基金（S2011040005355）

1.南方醫(yī)科大學(xué)公共衛(wèi)生與熱帶醫(yī)學(xué)學(xué)院生物統(tǒng)計(jì)學(xué)系（510515）

2.南方醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院數(shù)理系

△通信作者：歐春泉，E-mail：ouchunquan＠hotmail.com