三種模型對廣東省傷寒副傷寒逐月發病數預測的比較*

李 驪錢 俊楊 軍歐春泉△

三種模型對廣東省傷寒副傷寒逐月發病數預測的比較*

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

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

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

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

資料與方法

1.資料來源

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

2.統計方法及建模過程

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

（1）SARIMA模型［5］

SARIMA模型的構建分以下幾個步驟實現：

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

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

③參數估計，并檢驗參數是否有統計學意義。

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

（2）FARIMA模型

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

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

（3）DHR模型［6-7］

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

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

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

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

本研究DHR建模分以下幾個步驟實現：

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

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

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

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

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

（4）模型評價指標

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

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

（5）統計軟件

SARIMA和FARIMA建模采用R 2.15.1，DHR建模采用MATLAB 7.6編程實現。

結 果

1.發病趨勢分析

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

圖1 傷寒副傷寒發病數時序圖

2.各模型構建結果

（1）SARIMA模型

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

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

（2）FARIMA模型

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

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

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

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

（3）DHR模型

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

3.擬合預測效果評價

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

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

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

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

表2 三種模型的預測結果

討 論

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

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

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

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

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

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

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

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

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

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

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

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

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

16.胡世雄，邢慧嫻，鄧志紅.我國傳染病的預測預警現狀.中華預防醫學雜志，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.南方醫科大學公共衛生與熱帶醫學學院生物統計學系（510515）

2.南方醫科大學生物醫學工程學院數理系

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