葛申 馬建新 付凌姣 王晶 崔樹峰 張政

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

1 資料與方法

1.1資料來源 中國疾病預防控制系統(tǒng)2010年1月1日-2016年12月31日現住朝陽區(qū)的手足口病監(jiān)測數據。

1.2方法

1.2.1時間序列模型結構： 采用ARIMA模型，當觀測值為平穩(wěn)序列時，模型表達式為：

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

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

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

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

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

2 結果

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

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

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

圖2 原始時間序列經平穩(wěn)化后自相關系數函數圖(ACF圖)

圖3 原始時間序列經平穩(wěn)化后偏自相關系數函數圖(PACF圖)

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

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

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

圖5 2016年1-12月朝陽區(qū)手足口病月發(fā)病數趨勢預測圖

圖6 2016年7-12月朝陽區(qū)手足口病月發(fā)病數趨勢預測圖

3 討論

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

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