李 驪錢 俊楊 軍歐春泉△

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

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

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

SARIMA模型 FARIMA模型 DHR模型 預測 傷寒副傷寒

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

資料與方法

1.資料來源

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

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

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

（1）SARIMA模型［5］

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

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

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

③參數(shù)估計，并檢驗參數(shù)是否有統(tǒng)計學意義。

④模型診斷及預測，通過判斷殘差序列是否為白噪聲和存在自相關性評價模型是否合適。用確定的模型進行預測。

（2）FARIMA模型

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

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

（3）DHR模型［6-7］

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

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

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

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

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

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

②估計處理后序列的經驗譜。根據(jù)AIC最小化準則，自回歸譜階數(shù)設為24。

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

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

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

（4）模型評價指標

采用以下指標評價模型的預測準確性：均方根誤差（RMSE），平均絕對誤差（MAE），平均絕對百分比誤差（MAPE）及平均絕對偏差百分比（PMAD）。各指標的計算方法如下：

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

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

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

結 果

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

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

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

2.各模型構建結果

（1）SARIMA模型

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

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

（2）FARIMA模型

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

圖2 SARIMA（0，0，0）×（0，1，0）12殘差和SARIMA（1，0，0）×（0，1，0）12殘差的自相關和偏自相關圖

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

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

（3）DHR模型

模型殘差序列的ACF和PACF圖提示模型已經消除了序列的自相關特性，殘差序列通過了Ljung-Box檢驗，說明殘差為白噪聲。綜上，模型能較好地擬合該時間序列。

3.擬合預測效果評價

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

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

圖3 三種方法的預測結果圖

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

表2 三種模型的預測結果

討 論

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

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

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

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

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

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

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

5.彭志行，陶紅，賈成梅，等.時間序列分析在麻疹疫情預測預警中的應用研究.中國衛(wèi)生統(tǒng)計，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年江西省傷寒副傷寒流行病學分析.海峽預防醫(yī)學雜志，2009，15（2）：27-28.

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

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

11.吳偉慎，解曉華，單愛蘭，等.天津市1990-2005年傷寒副傷寒流行趨勢分析.現(xiàn)代預防醫(yī)學，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.春雅麗，黃中敏，錢文平.應用指數(shù)曲線預測傷寒副傷寒疫情.上海預防醫(yī)學雜志，2005，17（10）：466-467.

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

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

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

（責任編輯：劉 壯）

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

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

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

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

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