沈彭 魏峰

(安徽大學(xué)經(jīng)濟(jì)學(xué)院， 合肥 230039)

細(xì)菌性痢疾是新疆地區(qū)的高發(fā)傳染病，其發(fā)病率長(zhǎng)期位列該地區(qū)甲、乙類傳染病發(fā)病率前五名。防范細(xì)菌性痢疾仍是該地區(qū)疾病防控工作的重點(diǎn)之一。為了從根本上做好疾病防控工作，學(xué)者們開(kāi)始從傳染病的流行規(guī)律及預(yù)測(cè)方法等方面著手進(jìn)行研究。國(guó)內(nèi)已有學(xué)者利用時(shí)間序列模型對(duì)區(qū)域內(nèi)的疾病發(fā)生規(guī)律進(jìn)行預(yù)測(cè)，短期數(shù)據(jù)的預(yù)測(cè)效果較好[1-5]。本次研究擬利用新疆法定傳染病報(bào)告中細(xì)菌性痢疾的月發(fā)病率資料分析其時(shí)間變化趨勢(shì)，應(yīng)用ARIMA (Autoregressive Integrated Moving Average Model，即求和自回歸移動(dòng)平均模型)季節(jié)乘積模型對(duì)其進(jìn)行擬合和預(yù)測(cè)，探討將該模型用于細(xì)菌性痢疾疫情預(yù)測(cè)預(yù)警的可行性。

1 數(shù)據(jù)準(zhǔn)備與建模方法

1.1 數(shù)據(jù)準(zhǔn)備

數(shù)據(jù)為新疆地區(qū)2004 — 2016年細(xì)菌性痢疾月發(fā)病率，由中國(guó)疾病預(yù)防控制中心子網(wǎng)站 —— 公共衛(wèi)生科學(xué)數(shù)據(jù)中心披露[6]。其中，取2004年1月 — 2016年4月份的數(shù)據(jù)用于建立模型，取2016年4月至12月的數(shù)據(jù)用于預(yù)測(cè)值分析比較。ARIMA季節(jié)乘積模型對(duì)于數(shù)據(jù)的短期預(yù)測(cè)比較靈敏，所以只保留8個(gè)月的數(shù)據(jù)用于實(shí)際值與預(yù)測(cè)值的對(duì)比檢驗(yàn)。

1.2 ARIMA季節(jié)乘積模型的建立

ARIMA模型又稱Box-Jenkins模型，該分析方法是一種時(shí)間序列預(yù)測(cè)方法，用于描述非平穩(wěn)性資料。當(dāng)時(shí)間序列中存在各種趨勢(shì)性或季節(jié)波動(dòng)性時(shí)，可以建立連續(xù)模型ARIMA(p，d，q)或季節(jié)模型ARIMA(P，D，Q)S。有的季節(jié)性時(shí)間序列不僅含有季節(jié)性成分，還混有非季節(jié)性成分，此時(shí)可將兩者的混合效應(yīng)相乘，構(gòu)建ARIMA季節(jié)乘積模型ARIMA(p，d，q)(P，D，Q)S。其中，p、q分別為自回歸階數(shù)和移動(dòng)平均階數(shù)，P、Q分別為季節(jié)性自回歸階數(shù)和移動(dòng)平均階數(shù)，d為差分次數(shù)，D為季節(jié)性差分次數(shù)，S為季節(jié)周期步長(zhǎng)。模型結(jié)構(gòu)通常為：

(1)

式中：θ(B)為連續(xù)模型的q階移動(dòng)平均系數(shù)多項(xiàng)式；φ(B)為連續(xù)模型的p階自回歸系數(shù)多項(xiàng)式；θS(B)為季節(jié)模型的q階移動(dòng)平均系數(shù)多項(xiàng)式；φS(B)為季節(jié)模型的p階自回歸系數(shù)多項(xiàng)式。按以下流程建立模型(見(jiàn)圖1)[7]。

圖1 ARIMA模型建立流程

1.3 模型建立流程

(1) 數(shù)據(jù)預(yù)處理。在獲得一組觀察序列后，首先進(jìn)行平穩(wěn)性和白噪聲檢驗(yàn)。根據(jù)時(shí)序圖、(偏)自相關(guān)圖對(duì)數(shù)據(jù)的平穩(wěn)性作出初步判斷，并進(jìn)行PP單位根檢驗(yàn)(PP)、Dickey-Fuller單位根檢驗(yàn)(DF)或增廣Dickey-Fuller單位根檢驗(yàn)(ADF)。對(duì)于白噪聲檢驗(yàn)，則通過(guò)Q值或QLB統(tǒng)計(jì)量的P值進(jìn)行判斷。平穩(wěn)非白噪聲序列可進(jìn)入模型擬合階段，非平穩(wěn)序列需進(jìn)行平穩(wěn)化處理，白噪聲序列則不能進(jìn)行模型擬合。

(2) 數(shù)據(jù)平穩(wěn)化處理。通過(guò)一般差分與季節(jié)差分消除趨勢(shì)和季節(jié)的影響，并確定d、D值，差分后的序列需再次通過(guò)平穩(wěn)性和白噪聲的檢驗(yàn)。

(3) 模型識(shí)別。根據(jù)差分過(guò)程和結(jié)果，確定帶有季節(jié)因子的乘積模型ARIMA(p，d，q)(P，D，Q)S的基本形式，從低階到高階逐次完成模型的檢驗(yàn)識(shí)別。

(4) 參數(shù)估計(jì)與模型檢驗(yàn)。選用各參數(shù)有統(tǒng)計(jì)學(xué)意義以及AIC值、BIC值較小的適用模型進(jìn)行模型擬合，并對(duì)模型擬合結(jié)果進(jìn)行診斷，判斷殘差是否為白噪聲序列，從而確定模型的擬合優(yōu)度。

(5) 模型預(yù)測(cè)。利用完成擬合的模型進(jìn)行數(shù)據(jù)預(yù)測(cè)，并對(duì)預(yù)測(cè)結(jié)果進(jìn)行評(píng)價(jià)。

(6) 統(tǒng)計(jì)分析。應(yīng)用統(tǒng)計(jì)學(xué)軟件SAS 9.4完成數(shù)據(jù)處理、模型識(shí)別、模型預(yù)測(cè)等。

2 模型處理結(jié)果

2.1 數(shù)據(jù)預(yù)處理

圖2所示為新疆地區(qū)2004年1月 — 2016年4月細(xì)菌性痢疾月發(fā)病率時(shí)序圖。每年4月份的數(shù)據(jù)存在明顯的季節(jié)性變化趨勢(shì)，每年5月 — 9月發(fā)病率曲線處于波峰位置，為發(fā)病高峰期。圖3所示為2004年1月 — 2016年4月新疆細(xì)菌性痢疾月發(fā)病率自相關(guān)圖，其自相關(guān)系數(shù)也呈現(xiàn)出周期性的變化規(guī)律。

圖2 2004年1月 — 2016年4月新疆地區(qū)細(xì)菌性痢疾月發(fā)病率時(shí)序圖

2.2 數(shù)據(jù)平穩(wěn)化處理

根據(jù)數(shù)據(jù)特征，對(duì)原始序列進(jìn)行季節(jié)差分處理。差分后的自相關(guān)圖(見(jiàn)圖4)顯示，原序列的周期性趨勢(shì)消失； ADF檢驗(yàn)的結(jié)果(見(jiàn)圖5)顯示，差分過(guò)后P值均小于0.05，季節(jié)差分后的序列可以確定為平穩(wěn)序列。

2.3 白噪聲檢驗(yàn)

對(duì)季節(jié)差分后的序列延遲6、12、18、24期進(jìn)行白噪聲檢驗(yàn)(見(jiàn)圖6)，獲得的QLB統(tǒng)計(jì)量的P值均小于0.01，說(shuō)明差分后的序列是非隨機(jī)的。

圖3 2004年1月 — 2016年4月新疆細(xì)菌性痢疾月發(fā)病率自相關(guān)圖

圖4 2004年1月 — 2016年4月新疆細(xì)菌性痢疾月發(fā)病率季節(jié)差分后自相關(guān)圖

圖5 季節(jié)差分后ADF檢驗(yàn)結(jié)果

圖6 季節(jié)差分后白噪聲檢驗(yàn)的結(jié)果

2.4 模型識(shí)別

根據(jù)差分過(guò)程及其結(jié)果可知，季節(jié)差分1次，季節(jié)步長(zhǎng)為12，確定適合細(xì)菌性痢疾月發(fā)病率序列的模型初步形式為：ARIMA(p,0,q)×(P,1,Q)12。

首先對(duì)連續(xù)性模型的p、q進(jìn)行識(shí)別，由表1可知，BIC(0,1) = 0.814 058為最小。選擇ARIMA(0,0,1)為連續(xù)模型ARIMA(p,0,q)的最優(yōu)模型。

表1 不同p、q值情況下ARIMA(p,0,q)模型的BIC值

在連續(xù)模型的基礎(chǔ)上，對(duì)季節(jié)模型的P、Q進(jìn)行識(shí)別，從低階到高階逐個(gè)進(jìn)行嘗試，分別比較ARIMA(0,0,1)(0,1,0)12、ARIMA(0,0,1)(0,1,1)12、ARIMA(0,0,1)(1,1,0)12、ARIMA(0,0,1)(1,1,1)12、ARIMA(0,0,1)(0,1,2)12等模型的擬合優(yōu)度。在表2中，ARIMA(0,0,1)(0,1,0)12模型的AIC值和SBC值較小(AIC值為502.973 4,SBC值為508.798 7)，為相對(duì)最優(yōu)模型。

表2 不同P、Q值情況下ARIMA(P,1,Q)12模型的AIC值和SBC值

2.5 參數(shù)估計(jì)與模型檢驗(yàn)

如圖7所示，對(duì)常數(shù)項(xiàng)、MA(1，1)的估計(jì)值(-0.430 5、-0.653 4)進(jìn)行顯著性檢驗(yàn)，它們均具有統(tǒng)計(jì)學(xué)意義(P<0.05)，可以判定以下最終獲得的ARIMA(0,0,1)(0,1,0)12模型是恰當(dāng)?shù)摹?/p>

圖7 參數(shù)估計(jì)與模型檢驗(yàn)結(jié)果

2.6 模型預(yù)測(cè)

用ARIMA(0,0,1)(0,1,0)12模型預(yù)測(cè)新疆地區(qū)2016年5月 — 12月細(xì)菌性痢疾的月發(fā)病率，結(jié)果如圖8所示。圖中，擬合值與實(shí)際值基本重合，說(shuō)明模型對(duì)原始序列信息提取比較充分，擬合效果比較好；同時(shí)，預(yù)測(cè)值的變化趨勢(shì)與實(shí)際值基本一致，實(shí)際值基本落在預(yù)測(cè)曲線95%的可信區(qū)間內(nèi)。預(yù)測(cè)結(jié)果見(jiàn)表3。

圖8 ARIMA(0,0,1)(0,1,0)12模型擬合的新疆細(xì)菌性痢疾月發(fā)病率

表3 2016年5～12月新疆細(xì)菌性痢疾月發(fā)病率預(yù)測(cè)結(jié)果

3 討論

在實(shí)際應(yīng)用中，數(shù)據(jù)的時(shí)間序列大多是不平穩(wěn)的，表現(xiàn)出各種趨勢(shì)性或季節(jié)波動(dòng)性。有的季節(jié)性時(shí)間序列不僅含有季節(jié)性成分，還混有非季節(jié)性成分，若采用簡(jiǎn)單的時(shí)間序列模型進(jìn)行分析，往往不能得出準(zhǔn)確的預(yù)測(cè)結(jié)果。建立ARIMA季節(jié)乘積模型，可根據(jù)預(yù)測(cè)變量自身在各時(shí)刻變化規(guī)律的依存關(guān)系，以時(shí)間因素綜合替代各種影響因素，反映出疾病發(fā)生的趨勢(shì)性、周期性及季節(jié)性[8]。

本次研究中，以新疆地區(qū)2004年1月 — 2016年4月細(xì)菌性痢疾月報(bào)告發(fā)病率數(shù)據(jù)為基礎(chǔ)，擬合ARIMA季節(jié)乘積模型，同時(shí)對(duì)2016年5 — 12月的數(shù)據(jù)進(jìn)行了預(yù)測(cè)及比較。研究結(jié)果表明，采用ARIMA(0,0,1)(0,1,0)12模型擬合所得到的細(xì)菌性痢疾月發(fā)病率與實(shí)際報(bào)告數(shù)據(jù)基本重合，說(shuō)明模型對(duì)原始序列信息的提取比較充分；同時(shí)，預(yù)測(cè)值的變化趨勢(shì)與實(shí)際值變化一致，實(shí)際值基本落在預(yù)測(cè)曲線95%的可信區(qū)間內(nèi)。利用ARIMA季節(jié)乘積模型對(duì)新疆地區(qū)細(xì)菌性痢疾月報(bào)告發(fā)病率數(shù)據(jù)進(jìn)行擬合和預(yù)測(cè)，是合理可行的。

2010 — 2016年細(xì)菌性痢疾月報(bào)告發(fā)病率整體上低于2004 — 2009年的同期水平，特別是2012年最低，這可能與監(jiān)測(cè)方案的調(diào)整有關(guān)。因此，2012年發(fā)病高峰處預(yù)測(cè)值和真實(shí)值之間存在一定的差異，預(yù)測(cè)值要高于真實(shí)值。這一局限性決定了通過(guò)單次分析結(jié)果所建立的模型，需要不斷加入新的實(shí)際值，以修正或重新擬合才能更準(zhǔn)確地反映實(shí)際情況[9]。