胡建利 梁 祁 吳 瑩 劉文東 艾 靜 李 媛 張永杰 彭志行 鮑昌俊

傳染病的發(fā)病受到很多因素影響，而且影響因素之間又存在錯綜復雜的關系，因此靜態(tài)的因果結構模型很難揭示其流行趨勢，而根據事物自身變動規(guī)律建立動態(tài)模型——時間序列進行預測分析則是一種行之有效的方法。自回歸移動平均模型(Autoregressive Integrated Moving Average Model，ARIMA)是一類常用的隨機時間序列模型，已廣泛地應用于金融保險、社會科學、自然科學等領域中〔1，2〕，能比較準確地對序列未來各期進行預測。

腸道、呼吸道和自然疫源性傳染病的發(fā)病由于受到季節(jié)性因素或其他一些固有因素的影響，存在明顯的周期性變化。描述這類資料需要使用季節(jié)時間序列模型(Seasonal ARIMA Model，SARIMA)。本文運用SARIMA模型對江蘇省菌痢的月發(fā)病數建立數學模型，定量地預測其發(fā)病情況，并探討該模型進行傳染病預警的可行性。

資料和方法

1.資料來源

1990～2003年的菌痢統(tǒng)計數據來源于江蘇省法定傳染病年報表，2004～2010年的菌痢統(tǒng)計數據來源于江蘇省疾病監(jiān)測信息報告管理系統(tǒng)。

2.方法與原理

(1)基本思想

SARIMA模型:較早的文獻也稱其為乘積ARIMA模型，是隨機季節(jié)模型與ARIMA模型的結合，對于時間序列{Zt，t=1，2，…}有季節(jié)性、趨勢性和周期性時，可以建立非平穩(wěn)季節(jié)模型，表示為SARIMA(p，d，q)(P，D，Q)模型。其一般形式為〔3〕:

其中:φp(L)=1－φ1L－φ2L2－… －φpLp，p為非季節(jié)自回歸階數。

ΦP(Ls)=1－ΦsLs－Φ2sL2s－… －ΦPsLps，P為季節(jié)自回歸階數。

θq(L)=1－θ1L－θ2L2－… －θqLq，q為非季節(jié)移動平均階數。

ΘQ(Ls)=1－ΘsLs－Θ2sL2s－… －ΘQsLQs，Q 為季節(jié)自回歸階數。

d，D分別為普通差分和季節(jié)差分的階數，s為季節(jié)的長度。εt為白噪聲序列。

(2)建模過程

①數據的平穩(wěn)化

在確定時間序列模型之前需把不平穩(wěn)的時間序列轉化為平穩(wěn)的序列。通常將原序列進行自然對數變換消除其異方差，然后根據變換后序列的自相關和偏自相關圖，確定非季節(jié)差分階數d和季節(jié)差分階數D，d和D宜取較低階(通常取1，2，3)，s可以根據疾病的背景知識獲得。

②模型參數的估計

根據變換后的平穩(wěn)時間序列進行分析，尤其是序列的自相關和偏自相關圖，估計模型p、P、q、Q的值，采用最大似然估計或最小二乘法估計等對初步估計的模型進行檢驗。如果檢驗不通過，則調整參數，重新估計并檢驗，直至檢驗通過為止。

估計的模型通過檢驗是指:模型的參數必須通過t檢驗，且全部特征根的倒數都小于1〔4〕。

③模型的診斷檢驗

模型參數估計后，應該對模型的殘差是否為白噪聲進行檢驗，若殘差序列不是白噪聲序列，意味著殘差序列還存在有信息沒被提取，需要進一步改進模型。

實際運用中，可以獲得多個時序模型，為了得到一種最佳模型，可借助擬合優(yōu)度統(tǒng)計量來對比各個模型的優(yōu)劣。其中最常用的是調整后的決定系數、AIC和SC統(tǒng)計量。

3.數據分析

采用Eviews 5.0軟件進行數據的處理和分析。1990年1月～2009年12月菌痢月發(fā)病數數據用于建立模型，2010年1～5月數據用于驗證模型的預測效果。

結 果

1.資料分析

江蘇省1999年1月至2009年12月菌痢逐月發(fā)病數(Zt)曲線呈明顯的非平穩(wěn)性和季節(jié)性，并伴隨一定的周期性波動，見圖1。菌痢屬于腸道傳染病，發(fā)病有明顯的高峰季節(jié):每年7～8月份發(fā)病率最高，12月至次年2月發(fā)病率最低。

圖1 原始數據Zt序列圖

2.建立數學模型

(1)數據的平穩(wěn)化

從圖1可以看出，原始數據序列隨著時間呈現遞減型異方差。因此對原始數據首先進行自然對數轉換，以平穩(wěn)序列的方差。對數變換后的菌痢月發(fā)病數據(LnZt)自相關圖和偏自相關圖，見圖2。從圖2中可以看出自相關圖衰減很慢，說明LnZt是非平穩(wěn)的，且相關圖存在周期為12個月的季節(jié)波動。因此對LnZt進行一階非季節(jié)差分和一階季節(jié)差分，得到ΔΔ12LnZt。從序列ΔΔ12LnZt的相關圖和偏相關圖(圖3)可以看出，其自相關函數快速衰減，近似為一個平穩(wěn)過程。

(2)模型參數的估計

由于原始序列Zt對數變換后，經過一階非季節(jié)性差分和一階季節(jié)性差分達到平穩(wěn)，因此d=1、D=1。觀察序列ΔΔ12LnZt的偏相關圖，序列ΔΔ12LnZt的偏相關函數在滯后2階、12階、24階顯著地不為零(超過其95%的置信區(qū)間)，因此p=2(盡量選取低階);序列ΔΔ12LnZt的自相關函數在滯后2階、12階顯著地不為零，因此q=2;由于相關圖和偏自相關圖在滯后12階都顯著不為零，因此P=1、Q=1。

圖2 序列LnZt的相關圖(下)和偏相關圖(上)

圖3 序列ΔΔ12LnZt的相關圖(下)和偏相關圖(上)

首先考慮建立 SARIMA(2，1，2)(1，1，1)12模型，結果見表1，其中變量AR(2)的t=0.1917、P=0.8482＞0.05，MA(2)的 t= －1.0041、P=0.3165 ＞0.05，兩者都沒有通過t檢驗。然后刪除變量AR(2)，嘗試建立 SARIMA(1，1，2)(1，1，1)12模型，以及刪除變量MA(2)，嘗試建立 SARIMA(2，1，1)(1，1，1)12模型。表 1 可見，SARIMA(1，1，2)(1，1，1)12模型和 SARIMA(2，1，1)(1，1，1)12模型的所有參數都通過了 t檢驗。而且，SARIMA(1，1，2)(1，1，1)12模型和 SARIMA(2，1，1)(1，1，1)12模型均有 27 個根，包括 7 個實根和20個復根，其倒數均小于1。

(3)模型的診斷檢驗

對 SARIMA(2，1，1)(1，1，1)12模型殘差進行是否為白噪聲的Q統(tǒng)計量檢驗，該殘差序列的樣本量n為213，最大滯后期m可以取［213/10］或［213］，這里取22。Q22=18.534，P=0.356 ＞0.05，故不能拒絕殘差序列為白噪聲的原假設，檢驗通過。

對 SARIMA(1，1，2)(1，1，1)12模型殘差進行是否為白噪聲的Q統(tǒng)計量檢驗，該殘差序列的樣本量n為214，最大滯后期m可以?。?14/10］或［214］，這里取22。Q22=17.662，P=0.410 ＞0.05，故不能拒絕殘差序列為白噪聲的原假設，檢驗通過。

從表 2 可以看出，SARIMA(1，1，2)(1，1，1)12模型的 R2和調整 R2均比 SARIMA(2，1，1)(1，1，1)12的大，且對于AIC統(tǒng)計量和SC統(tǒng)計量，SARIMA(1，1，2)(1，1，1)12模型都比 SARIMA(2，1，1)(1，1，1)12模型要小。因此，SARIMA(1，1，2)(1，1，1)12模型擬合效果較好。

表1 各種SARIMA模型的檢驗結果

表2 兩種模型的擬合優(yōu)度比較

3.模型的預測

SARIMA(1，1，2)(1，1，1)12模型的表達式為:

(1－0.6249L)(1－0.2113L12)(1－L)(1－L12)log(Zt)=(1－0.6899L－0.2355L2)(1－0.9198L12)εt。該模型預測擬合圖(見圖4)顯示，實際數據與預測數據相當吻合;對1990年1月至2009年12月的菌痢發(fā)病數進行回代預測，結果顯示平均誤差率為13.89%。

圖 4 SARIMA(1，1，2)(1，1，1)12模型預測擬合圖

根據所建模型對2010年1～5月的菌痢發(fā)病數進行短期預測，預測結果分別為:304、217、329、390和598，此5 個月的實際發(fā)病數分別為:277、262、268、414和601，其預測誤差率分別為:9.75%、－17.18%、22.76%、－5.80%和－0.50%。

討 論

目前，對傳染病發(fā)病率進行預測時，常用的模型有曲線擬合、灰色模型、Markov模型、ARIMA/SARIMA模型等。其中，ARIMA/SARIMA模型將傳染病流行過程中各種影響因素的綜合效應統(tǒng)一蘊涵于時間序列中，這是其應用于傳染病預測的一個突出優(yōu)點。腸道、呼吸道和自然疫源性傳染病的發(fā)病情況有著明顯的季節(jié)性和周期性，如果不考慮這些因素的影響，做出的預測往往不準確，因此SARIMA模型在該領域有著廣泛的適用性。該模型的建立已經有一套明確的準則，適用于各種復雜的時序模式，同時一些統(tǒng)計軟件(如SAS、SPSS、Eviews等)為該模型的建立提供了有利條件，使其得到了廣泛的應用。本文利用SARIMA模型對江蘇省菌痢未來疫情動態(tài)和發(fā)展趨勢進行預測，得到較好的預測效果，對存在季節(jié)性、周期性波動的傳染病的預警具有指導意義。

江蘇省菌痢的發(fā)病率從20世紀90代初的100.00/10萬以上，下降至2009年的9.96/10萬，但它仍居腸道傳染病發(fā)病率的首位。通過SARIMA模型對江蘇省菌痢2010年1～5月的月發(fā)病數進行了預測，發(fā)現我省菌痢的發(fā)病將繼續(xù)呈下降趨勢。但由于流行因素的廣泛存在，作為一種高發(fā)病率的疾病，菌痢的發(fā)病率仍將維持在一個較高的水平，應繼續(xù)實施現行防制策略，加強病原學監(jiān)測、疾病預警工作。

ARIMA模型或SARIMA模型對疾病進行預測分析，有兩大優(yōu)點:①利用預測變量自身的變化規(guī)律建立模型，不考慮其相關因素;②明確考慮時間序列的非平穩(wěn)性，通過取對數、差分等方法把序列平穩(wěn)化后，再考慮建模問題。但是應用該模型應該注意:①至少需要50個以上的歷史數據;②所建立的模型，不能作為永久不變的預測工具，只能用于短期預測。對于已建立的模型應不斷加入新的實際值，以修正或重新擬合更優(yōu)的模型〔6〕。

1．Helfenstein U．Box-Jenkins modelling in medical research．Stat Methods Med Res，1996，5(1):3-22．

2．Kao JJ，Huang SS．Forecasts using neural network versus Box-Jenkins methodology for ambient air quality monitoring data．J Air Waste Manag Assoc，2000，50(2):219-226．

3．樊歡歡，張凌云主編．Eviews統(tǒng)計分析與應用．北京:機械工作出版社，2009:227-228．

4．中國人民銀行調查統(tǒng)計司主編．時間序列 X-12-ARIMA季節(jié)調整——原理與方法．北京:中國金融出版社，2006:59-60．

5．易丹輝主編．數據分析與EVIEWS應用．北京:中國人民大學出版社，2008:137-140．

6．漆莉，李革，李勤．ARIMA模型在流行性感冒預測中的應用．第三軍醫(yī)大學學報，2007，29(3):267-269．