葛申 馬建新 付凌姣 王晶 崔樹(shù)峰 張政

手足口病是由多種腸道病毒引起的,以發(fā)熱及手、足、口、臀等部位特征性皮疹或皰疹為主要癥狀的兒童常見(jiàn)傳染病[1]。2008年，安徽省阜陽(yáng)市曾出現(xiàn)較大規(guī)模流行[2]，同年5月，我國(guó)將其列入法定報(bào)告?zhèn)魅静」芾怼＝鼛啄?，我?guó)手足口病呈現(xiàn)發(fā)病強(qiáng)度高、高峰持續(xù)時(shí)間長(zhǎng)、疫情分布廣等特點(diǎn)[3]。北京市屬于“熱點(diǎn)”發(fā)病城市[4]，朝陽(yáng)區(qū)手足口病發(fā)病數(shù)一直位于北京市發(fā)病數(shù)排名前列[5-6]。手足口病發(fā)病受到許多復(fù)雜因素的影響，一般回歸預(yù)測(cè)很難獲取并分析全部的相關(guān)因素。時(shí)間序列分析利用事物發(fā)展延續(xù)性規(guī)律，以時(shí)間為自變量，用歷史監(jiān)測(cè)值建立模型預(yù)測(cè)事物未來(lái)情況[7]。自回歸移動(dòng)平均模型(autoregressive integrated moving average model，ARIMA)是重要的時(shí)間序列分析預(yù)測(cè)模型，其操作方法方便、獲取數(shù)據(jù)成本低、適用范圍廣、實(shí)用性強(qiáng)、短期預(yù)測(cè)精度較高，在傳染病預(yù)測(cè)、預(yù)警方面有較多的應(yīng)用。在預(yù)測(cè)手足口病方面，有研究表明ARIMA模型可以較好預(yù)測(cè)短期內(nèi)的變化趨勢(shì)[8]。因此，現(xiàn)通過(guò)ARIMA模型對(duì)朝陽(yáng)區(qū)手足口病的發(fā)病進(jìn)行時(shí)間序列分析并建立預(yù)測(cè)模型，以研究朝陽(yáng)區(qū)手足口病流行特征和發(fā)展規(guī)律，為朝陽(yáng)區(qū)手足口病的防控工作提供科學(xué)的參考依據(jù)。

1 資料與方法

1.1資料來(lái)源 中國(guó)疾病預(yù)防控制系統(tǒng)2010年1月1日-2016年12月31日現(xiàn)住朝陽(yáng)區(qū)的手足口病監(jiān)測(cè)數(shù)據(jù)。

1.2方法

1.2.1時(shí)間序列模型結(jié)構(gòu)： 采用ARIMA模型，當(dāng)觀測(cè)值為平穩(wěn)序列時(shí)，模型表達(dá)式為：

Yt=φ1Yt-1+φ2Yt-2+…+φpYt-p+

et-θ1et-1-θ2et-2-…-θqet-q

模型中Yt是時(shí)間序列在t期的觀測(cè)值。當(dāng)序列不平穩(wěn)時(shí)，通過(guò)差分使其平穩(wěn)，模型為ARIMA (p,d,q)。模型參數(shù)：p為自回歸階數(shù)，d為差分次數(shù)、q為移動(dòng)平均階數(shù)。若進(jìn)一步考慮資料的季節(jié)性/周期性，則模型標(biāo)記為ARIMA (p,d,q)(P,D,Q)s。新的參數(shù)P、Q為季節(jié)自回歸和移動(dòng)平均的階數(shù)，D為季節(jié)性差分次數(shù)，s為季節(jié)性周期循環(huán)的長(zhǎng)度，本文通過(guò)月發(fā)病數(shù)預(yù)測(cè)以12個(gè)月為季節(jié)性循環(huán)周期。

1.2.2模型建模步驟： 主要包括四個(gè)階段：①序列的平穩(wěn)化：ARIMA模型應(yīng)用條件是預(yù)測(cè)對(duì)象的時(shí)間序列為平穩(wěn)隨機(jī)序列，不平穩(wěn)的序列需要進(jìn)行預(yù)處理，處理后通過(guò)進(jìn)行分析確認(rèn);②序列特征識(shí)別和模型的識(shí)別：通過(guò)繪制并觀察繪制時(shí)間序列圖、自相關(guān)系數(shù)函數(shù)圖(ACF圖)和偏自相關(guān)系數(shù)函數(shù)圖(PACF圖)進(jìn)行模型的初步識(shí)別，對(duì)模型進(jìn)行定階;③參數(shù)估計(jì)和模型診斷：利用非線性最小二乘法估計(jì)模型參數(shù)，采取從低階到高階逐步嘗試的方法，依次擬合不同參數(shù)組合。確定參數(shù)后，對(duì)原始數(shù)據(jù)與擬合數(shù)據(jù)的殘差序列進(jìn)行白噪聲檢測(cè)，其檢驗(yàn)方法為計(jì)算Box-Jenkins統(tǒng)計(jì)量(Q值);④模型預(yù)測(cè)應(yīng)用：選定最佳模型后，對(duì)2016年朝陽(yáng)區(qū)月發(fā)病數(shù)進(jìn)行預(yù)測(cè)，將預(yù)測(cè)發(fā)病數(shù)與實(shí)際發(fā)病數(shù)進(jìn)行比較，以驗(yàn)證模型效果。

1.3統(tǒng)計(jì)學(xué)方法 使用Excel 2013軟件對(duì)收集到的數(shù)據(jù)進(jìn)行匯總、整理，使用SPSS 21.0軟件對(duì)資料進(jìn)行統(tǒng)計(jì)分析和模型構(gòu)建。

2 結(jié)果

2.12010-2015年朝陽(yáng)區(qū)手足口病發(fā)病時(shí)間趨勢(shì) 通過(guò)繪制2010-2015年手足口病月發(fā)病數(shù)時(shí)間序列圖(圖1)進(jìn)行分析：朝陽(yáng)區(qū)手足口病各月均有發(fā)病，呈明顯季節(jié)性規(guī)律，5-7月出現(xiàn)明顯的流行高峰，1-2月發(fā)病數(shù)最低，10-11月出現(xiàn)發(fā)病數(shù)反彈的情況，特別是2011年和2014年出現(xiàn)較為明顯的小高峰，其他年份此特征不顯著。從流行長(zhǎng)期趨勢(shì)來(lái)看，高發(fā)年份之后次年發(fā)病大幅下降，此后逐年小幅升高。

圖1 2010-2015年朝陽(yáng)區(qū)手足口病發(fā)病數(shù)時(shí)序分布

2.2模型構(gòu)建 通過(guò)前述發(fā)病趨勢(shì)分析發(fā)現(xiàn)，2010-2015年北京市朝陽(yáng)區(qū)手足口病發(fā)病時(shí)間序列存在以12個(gè)月為1個(gè)周期的季節(jié)性特征，且序列波動(dòng)較大，是不平穩(wěn)的序列。對(duì)序列進(jìn)行季節(jié)性差分和自然對(duì)數(shù)轉(zhuǎn)換以平穩(wěn)季節(jié)性波動(dòng)和減少方差波動(dòng)。平穩(wěn)化后，通過(guò)觀察ACF圖(圖2)和PACF圖(圖3)，ACF圖呈震蕩衰減形式，PACF圖第一階函數(shù)值特征顯示，是典型的自回歸過(guò)程，確定模型為季節(jié)乘積模型ARIMA(p,d,q)(P,D,Q)12。

圖2 原始時(shí)間序列經(jīng)平穩(wěn)化后自相關(guān)系數(shù)函數(shù)圖(ACF圖)

圖3 原始時(shí)間序列經(jīng)平穩(wěn)化后偏自相關(guān)系數(shù)函數(shù)圖(PACF圖)

2.3參數(shù)估計(jì)和模型診斷 對(duì)模型的參數(shù)采取從低到高逐個(gè)進(jìn)行嘗試的辦法，得出不同階數(shù)組合的模型。經(jīng)過(guò)篩選后選取標(biāo)準(zhǔn)化BIC值最小，R2較大且簡(jiǎn)潔的模型為最佳模型。模型ARIMA(1,0,0)(1,1,0)12的各參數(shù)值均有意義，R2=0.679表明模型擬合程度較好，標(biāo)準(zhǔn)化BIC值=10.894，在備選的模型中最小。采用Ljung-Box方法做殘差白噪聲檢測(cè)，其Ljung-Box Q=22.59，P>0.05。做殘差序列的ACF、PACF圖(圖4)，殘差序列的自相關(guān)系數(shù)和偏自相關(guān)系數(shù)均落入95%CI中。證明殘差序列為隨機(jī)誤差，提示模型已經(jīng)將時(shí)間序列中蘊(yùn)含的信息提取出來(lái)。

圖4 模型ARIMA(1,0,0)(1,1,0)12殘差序列ACF、PACF圖

2.4模型預(yù)測(cè)應(yīng)用 通過(guò)前述建模方法，對(duì)朝陽(yáng)區(qū)2010年1月-2015年12月手足口病月發(fā)病數(shù)建模，對(duì)2016年1-12月的月發(fā)病數(shù)進(jìn)行預(yù)測(cè)(圖5)。以2016年1-12月的實(shí)際月發(fā)病數(shù)驗(yàn)證：預(yù)測(cè)結(jié)果與實(shí)際情況總體趨勢(shì)一致，季節(jié)性規(guī)律基本相同，實(shí)際發(fā)病數(shù)在預(yù)測(cè)發(fā)病數(shù)的95%CI范圍內(nèi)波動(dòng)，但預(yù)測(cè)發(fā)病數(shù)的平均相對(duì)誤差達(dá)49.37%。將2016年1-6月的發(fā)病數(shù)序列繼續(xù)納入模型，對(duì)2016年7-12月的發(fā)病數(shù)進(jìn)行預(yù)測(cè)(圖6)。再以2016年7-12月實(shí)際發(fā)病數(shù)為驗(yàn)證數(shù)據(jù)。結(jié)果可見(jiàn)，預(yù)測(cè)結(jié)果與實(shí)際情況的總體趨勢(shì)一致，平均相對(duì)誤差降低至18.12%，較之預(yù)測(cè)全年發(fā)病數(shù)相對(duì)誤差明顯減小，與ARIMA模型短期預(yù)測(cè)效果更好的特點(diǎn)相符。

圖5 2016年1-12月朝陽(yáng)區(qū)手足口病月發(fā)病數(shù)趨勢(shì)預(yù)測(cè)圖

圖6 2016年7-12月朝陽(yáng)區(qū)手足口病月發(fā)病數(shù)趨勢(shì)預(yù)測(cè)圖

3 討論

3.1此次擬合的ARIMA模型，預(yù)測(cè)手足口月發(fā)病趨勢(shì)與實(shí)際發(fā)病整體變化趨勢(shì)一致。適用于手足口病等具有季節(jié)性變動(dòng)特征的傳染病預(yù)測(cè)[9-11]。ARIMA模型通常用來(lái)處理平穩(wěn)的時(shí)間序列，而傳染病發(fā)病數(shù)據(jù)序列大多是非平穩(wěn)的。因此，建模之前需對(duì)數(shù)據(jù)進(jìn)行處理，以達(dá)到平穩(wěn)化的要求。分析發(fā)現(xiàn)，朝陽(yáng)區(qū)手足口病發(fā)病時(shí)間序列呈現(xiàn)出明顯的季節(jié)性特征，且每個(gè)月之間發(fā)病數(shù)相差較大，是非平穩(wěn)序列。因此，采用自然對(duì)數(shù)轉(zhuǎn)換和季節(jié)性差分減少方差波動(dòng)和實(shí)現(xiàn)序列平穩(wěn)化。由于時(shí)間序列是根據(jù)事物發(fā)展的延續(xù)性而建立的，克服了影響預(yù)測(cè)對(duì)象的因素錯(cuò)綜復(fù)雜、不易分析、數(shù)據(jù)資料不易獲取等問(wèn)題，以時(shí)間(t)綜合替代各種影響因素，其模型構(gòu)建短期精確度較高[12]。分析結(jié)果顯示，對(duì)2016年全年(1-12月)和下半年(7-12月)的發(fā)病數(shù)據(jù)進(jìn)行預(yù)測(cè)，實(shí)際發(fā)病數(shù)與預(yù)測(cè)發(fā)病數(shù)的平均相對(duì)誤差分別為49.37%和18.12%。因此，在實(shí)際應(yīng)用中，應(yīng)動(dòng)態(tài)開(kāi)展短期預(yù)測(cè)，將新的實(shí)際發(fā)病數(shù)據(jù)納入模型分析，從而提高預(yù)測(cè)精度，這也符合ARIMA模型依賴事物發(fā)展延續(xù)性的基本思想。

3.2ARIMA模型是依照事物發(fā)展的慣性趨勢(shì)預(yù)測(cè)未來(lái)發(fā)病趨勢(shì)的一種時(shí)間序列分析模型，而在實(shí)際中，影響手足口病發(fā)病的因素復(fù)雜且不斷變化，如：易感者的數(shù)量、氣候、病原譜的改變、手足口病防控工作的開(kāi)展、疫苗的使用等均會(huì)對(duì)疾病的發(fā)展趨勢(shì)造成影響[13-15]。但實(shí)際防病工作不能簡(jiǎn)單地依靠模型來(lái)判斷，需要結(jié)合流行病學(xué)專業(yè)理論知識(shí)及發(fā)病影響因素進(jìn)行具體分析，才能使模型預(yù)測(cè)在防控工作中發(fā)揮更大的作用。