季節(jié)時(shí)間序列模型在菌痢發(fā)病預(yù)測(cè)中的應(yīng)用

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

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

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

資料和方法

1.資料來(lái)源

1990～2003年的菌痢統(tǒng)計(jì)數(shù)據(jù)來(lái)源于江蘇省法定傳染病年報(bào)表，2004～2010年的菌痢統(tǒng)計(jì)數(shù)據(jù)來(lái)源于江蘇省疾病監(jiān)測(cè)信息報(bào)告管理系統(tǒng)。

2.方法與原理

(1)基本思想

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

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

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

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

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

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

(2)建模過(guò)程

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

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

②模型參數(shù)的估計(jì)

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

估計(jì)的模型通過(guò)檢驗(yàn)是指:模型的參數(shù)必須通過(guò)t檢驗(yàn)，且全部特征根的倒數(shù)都小于1〔4〕。

③模型的診斷檢驗(yàn)

模型參數(shù)估計(jì)后，應(yīng)該對(duì)模型的殘差是否為白噪聲進(jìn)行檢驗(yàn)，若殘差序列不是白噪聲序列，意味著殘差序列還存在有信息沒(méi)被提取，需要進(jìn)一步改進(jìn)模型。

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

3.數(shù)據(jù)分析

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

結(jié) 果

1.資料分析

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

圖1 原始數(shù)據(jù)Zt序列圖

2.建立數(shù)學(xué)模型

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

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

(2)模型參數(shù)的估計(jì)

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

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

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

首先考慮建立 SARIMA(2，1，2)(1，1，1)12模型，結(jié)果見表1，其中變量AR(2)的t=0.1917、P=0.8482＞0.05，MA(2)的 t= －1.0041、P=0.3165 ＞0.05，兩者都沒(méi)有通過(guò)t檢驗(yàn)。然后刪除變量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模型的所有參數(shù)都通過(guò)了 t檢驗(yàn)。而且，SARIMA(1，1，2)(1，1，1)12模型和 SARIMA(2，1，1)(1，1，1)12模型均有 27 個(gè)根，包括 7 個(gè)實(shí)根和20個(gè)復(fù)根，其倒數(shù)均小于1。

(3)模型的診斷檢驗(yàn)

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

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

從表 2 可以看出，SARIMA(1，1，2)(1，1，1)12模型的 R2和調(diào)整 R2均比 SARIMA(2，1，1)(1，1，1)12的大，且對(duì)于AIC統(tǒng)計(jì)量和SC統(tǒng)計(jì)量，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模型的檢驗(yàn)結(jié)果

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

3.模型的預(yù)測(cè)

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

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

圖 4 SARIMA(1，1，2)(1，1，1)12模型預(yù)測(cè)擬合圖

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

討 論

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

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

ARIMA模型或SARIMA模型對(duì)疾病進(jìn)行預(yù)測(cè)分析，有兩大優(yōu)點(diǎn):①利用預(yù)測(cè)變量自身的變化規(guī)律建立模型，不考慮其相關(guān)因素;②明確考慮時(shí)間序列的非平穩(wěn)性，通過(guò)取對(duì)數(shù)、差分等方法把序列平穩(wěn)化后，再考慮建模問(wèn)題。但是應(yīng)用該模型應(yīng)該注意:①至少需要50個(gè)以上的歷史數(shù)據(jù);②所建立的模型，不能作為永久不變的預(yù)測(cè)工具，只能用于短期預(yù)測(cè)。對(duì)于已建立的模型應(yīng)不斷加入新的實(shí)際值，以修正或重新擬合更優(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)計(jì)分析與應(yīng)用．北京:機(jī)械工作出版社，2009:227-228．

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

5．易丹輝主編．?dāng)?shù)據(jù)分析與EVIEWS應(yīng)用．北京:中國(guó)人民大學(xué)出版社，2008:137-140．

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