自回歸移動平均模型在北京市朝陽區手足口病發病預測中的應用

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

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

1 資料與方法

1.1資料來源 中國疾病預防控制系統2010年1月1日-2016年12月31日現住朝陽區的手足口病監測數據。

1.2方法

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

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

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

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

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

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

2 結果

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

圖1 2010-2015年朝陽區手足口病發病數時序分布

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

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

圖3 原始時間序列經平穩化后偏自相關系數函數圖(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模型預測應用 通過前述建模方法，對朝陽區2010年1月-2015年12月手足口病月發病數建模，對2016年1-12月的月發病數進行預測(圖5)。以2016年1-12月的實際月發病數驗證：預測結果與實際情況總體趨勢一致，季節性規律基本相同，實際發病數在預測發病數的95%CI范圍內波動，但預測發病數的平均相對誤差達49.37%。將2016年1-6月的發病數序列繼續納入模型，對2016年7-12月的發病數進行預測(圖6)。再以2016年7-12月實際發病數為驗證數據。結果可見，預測結果與實際情況的總體趨勢一致，平均相對誤差降低至18.12%，較之預測全年發病數相對誤差明顯減小，與ARIMA模型短期預測效果更好的特點相符。

圖5 2016年1-12月朝陽區手足口病月發病數趨勢預測圖

圖6 2016年7-12月朝陽區手足口病月發病數趨勢預測圖

3 討論

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

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