基于小波分析的手足口病時間序列預(yù)測*

余功超 馮慧芬 封 爽 趙 敬 徐 晶

【提 要】 目的 構(gòu)建基于小波分析的自回歸移動平均(ARIMA)模型預(yù)測手足口病流行，提高預(yù)測精度。方法 使用2010-2015年鄭州市疾控中心手足口病監(jiān)測數(shù)據(jù)，構(gòu)建基于小波分析的ARIMA模型進(jìn)行預(yù)測，用2016年數(shù)據(jù)進(jìn)行驗證，并與單純的ARIMA模型進(jìn)行比較。結(jié)果 構(gòu)建的基于小波分解一層的ARIMA模型為ARIMA(0，1，3)(2，1，0)52，矯正后的AIC=2747.82，殘差序列的ACF、PACF圖示殘差序列無自相關(guān)，Box-Ljung test統(tǒng)計量為0.9177，P=0.34，認(rèn)為該殘差為白噪聲序列，模型擬合良好。預(yù)測2016年發(fā)病趨勢與實際較為相符，均方根誤差RMSE(root mean square error)、平均絕對誤差MAE(mean absolute error)、平均絕對百分比誤差MAPE(mean absolute percentage error )，分別為49.42、26.45、15.75(訓(xùn)練集擬合)和275.84、219.90、72.95(驗證集預(yù)測)，除了驗證集MAPE外，均小于單一的ARIMA模型。結(jié)論 基于小波分析的ARIMA模型可用于手足口病時間序列預(yù)測，擬合和預(yù)測性能較單一的ARIMA模型好。

手足口病是一種由多種腸道病毒引起的傳染病[1]，好發(fā)于5歲以下兒童，是兒童主要感染性疾病之一。EV71和CA16是引起手足口病主要的病原，2015年EV71疫苗問世，有報道稱近年來EV71和CA16感染病例數(shù)有所下降，但其他腸道病毒感染病例數(shù)上升[2-3]，其發(fā)病率仍居高不下[4-5]，給社會造成了很大的醫(yī)療和經(jīng)濟(jì)負(fù)擔(dān)。手足口病無特殊治療方法，主要是對癥治療和防治并發(fā)癥[6]。因此，明晰手足口病流行周期及模式，更精準(zhǔn)地預(yù)測其流行，對衛(wèi)生行政部門制定預(yù)防控制策略、降低其發(fā)病率、減輕疾病負(fù)擔(dān)具有重要意義。

自回歸移動平均模型(ARIMA模型)是最常用的手足口病時間序列預(yù)測模型，適用于穩(wěn)定性時間序列，但現(xiàn)實世界的時間序列往往是非穩(wěn)定性的[7]，特別是流行病學(xué)時間序列[8]。雖然可通過各種變換和差分將之變成穩(wěn)定性時間序列，進(jìn)而達(dá)到應(yīng)用ARIMA模型的條件，但這種單一函數(shù)變換能力有限，需要技巧，仍然可能損失了一些有用的非線性信息[9]。小波分析是一種有效的適用于非穩(wěn)定性時間序列的時頻分析方法[8]。一些學(xué)者將小波分析和ARIMA模型結(jié)合起來進(jìn)行預(yù)測，研究結(jié)果表明結(jié)合了小波分析的ARIMA模型預(yù)測精度要比單一的ARIMA模型高。這種方法已經(jīng)應(yīng)用到建筑物沉降預(yù)報[10]，也有應(yīng)用到傳染病時間序列預(yù)測[9，11]。

基于此，本研究使用鄭州市疾控中心監(jiān)測數(shù)據(jù)，構(gòu)建基于小波分析的ARIMA模型進(jìn)行預(yù)測，以期提高手足口病時間序列預(yù)測精度，為手足口病防控提供參考。

資料與方法

1.數(shù)據(jù)來源

2010年至2016年鄭州市手足口病發(fā)病例數(shù)資料來自鄭州市疾控中心。將數(shù)據(jù)分為兩部分，2010年至2015年數(shù)據(jù)作為訓(xùn)練集，用來建模；2016年數(shù)據(jù)作為驗證集，用來檢驗?zāi)Ｐ皖A(yù)測性能。

2.ARIMA模型構(gòu)建

自回歸移動平均模型如果包含季節(jié)性則記為SARIMA(p，d，q)(P，D，Q)n，其中p，q為非季節(jié)性模型的自回歸及移動平均參數(shù)，d為普通差分的階數(shù)；P，Q為季節(jié)性模型的自回歸及移動平均參數(shù)，D為季節(jié)差分的階數(shù)，n為周期長度。建模步驟：(1)畫出時間序列圖，觀察趨勢變化，是否為平穩(wěn)性序列，是否有季節(jié)性周期性趨勢；(2)對非平穩(wěn)性時間序列進(jìn)行普通差分和季節(jié)性差分，將其變?yōu)榉€(wěn)定性時間序列，根據(jù)自相關(guān)系數(shù)、偏自相關(guān)系數(shù)及augmented Dickey-Fuller test結(jié)果(檢驗水準(zhǔn)α=0.05)判斷差分后的時間序列是否為穩(wěn)定性序列；(3)主要根據(jù)自相關(guān)系數(shù)和偏自相關(guān)系數(shù)，結(jié)合經(jīng)驗不斷嘗試確定模型的自回歸、移動平均參數(shù)。本研究使用R軟件“forecast”包中的自動建模功能(“auto.arima”)，可自動嘗試不同的自回歸、移動平均參數(shù)值來擬合模型，并選擇矯正的AIC(akaike′ information criterion)值最小的模型作為最優(yōu)模型(更多自動建模規(guī)則詳見Rob J Hyndman和George Athanasopoulos所著的《Forecasting：Principles and Practice》在線圖書，http：//otexts.com/fpp2/)；(4)作出模型殘差的自相關(guān)系數(shù)ACF圖、偏自相關(guān)系數(shù)PACF圖，判斷殘差是否具有自相關(guān)。使用Box-Ljung test 檢驗殘差是否為白噪聲序列，檢驗水準(zhǔn)α=0.05。若擬合效果好，模型提取信息充分，則殘差應(yīng)為白噪聲序列，無自相關(guān)，即該模型可用于預(yù)測。建模過程在R 3.6.3軟件中完成。

3.基于小波分析的ARIMA模型構(gòu)建

本研究使用離散小波分解，選擇的小波為Daubechies小波。設(shè)原始序列為X，利用Daubechies小波經(jīng)j層離散小波分解將原始序列分解為近似成分(approximation component)和細(xì)節(jié)成分(detail components)，使用小波重構(gòu)得到近似成分序列cA和細(xì)節(jié)成分序列cDj。則有

X=cA+cD1+cD2+…+cDj

其中近似成分序列是低頻的，與原始序列的輪廓較為一致，但比之平滑，細(xì)節(jié)成分序列是高頻的，通常是含有噪聲的細(xì)小的波動。我們分別用小波對原始序列進(jìn)行一層和兩層分解，使用MATLAB 軟件(Version R2014a)完成。用重構(gòu)的近似成分序列參照SARIMA模型構(gòu)建時最優(yōu)模型參數(shù)來構(gòu)建SARIMA模型。

4.模型評價標(biāo)準(zhǔn)

采用均方根誤差RMSE(root mean square error)、平均絕對誤差MAE(mean absolute error)、平均絕對百分比誤差MAPE(mean absolute percentage error )來評估模型擬合及預(yù)測性能。指標(biāo)的計算公式如下：

結(jié) 果

1.手足口病周發(fā)病例數(shù)時間序列圖

從圖1中可以看出手足口病周發(fā)病例數(shù)時間序列不平穩(wěn)，呈現(xiàn)比較明顯的季節(jié)性周期性特征，周期為1年(52周)，每年5月至7月達(dá)到高峰，部分年份還有冬季的小高峰。

圖1 2010-2016年鄭州市手足口病周發(fā)病例數(shù)的時間序列圖

2.SARIMA模型

對2010-2015年鄭州市手足口病發(fā)病例數(shù)時間序列進(jìn)行一階季節(jié)性差分(周期長度為52周)和一階差分后，其ACF圖和PACF圖示，序列已變?yōu)槠椒€(wěn)性序列(圖2)。ADF單位根檢驗結(jié)果：Dickey-Fuller=-6.9746，P=0.01，認(rèn)為差分后的序列為穩(wěn)定性序列。使用自動建模功能，我們得到最優(yōu)模型為SARIMA(0，1，3)(2，1，0)52，矯正AIC值(AICc)=2967.27。殘差序列的ACF、PACF圖示殘差無自相關(guān)(圖3)，Box-Ljung 檢驗統(tǒng)計量為0.0025，P=0.96，認(rèn)為該殘差為白噪聲序列，該模型擬合良好，可以用于預(yù)測。

圖2 一階季節(jié)性差分和普通差分后序列的ACF、PACF圖

圖3 SARIMA模型殘差序列ACF、PACF圖

3.基于小波分析的SARIMA模型

用小波分析對2010-2015年鄭州市手足口病發(fā)病例數(shù)時間序列進(jìn)行一層分解后，構(gòu)建的SARIMA模型為SARIMA(0，1，3)(2，1，0)52，AICc=2747.82，殘差序列的ACF、PACF圖示殘差序列無自相關(guān)(圖4)，Box-Ljung檢驗統(tǒng)計量為0.9177，P=0.34，認(rèn)為該殘差序列為白噪聲序列。用小波分析對原序列進(jìn)行兩層分解后，構(gòu)建的SARIMA模型為SARIMA(0，1，3)(2，1，0)52，AICc=2726.37，殘差序列的ACF、PACF圖示殘差序列無自相關(guān)(圖5)，Box-Ljung test 統(tǒng)計量為3.5971，P=0.06，認(rèn)為該殘差序列為白噪聲序列。兩模型均擬合良好，可以用于預(yù)測。

圖4 基于小波的SARIMA模型(分解一層)殘差序列ACF、PACF圖

圖5 基于小波的SARIMA模型(分解兩層)殘差序列ACF、PACF圖

6.模型比較

SARIMA模型、基于小波分解一層或兩層的SARIMA模型的擬合和預(yù)測圖見圖6，各模型擬合和預(yù)測與實際都較為相符。三種模型的擬合與預(yù)測性能指標(biāo)見表1，無論是擬合還是預(yù)測，基于小波分析的SARIMA模型(分解一層)的各項指標(biāo)均小于小波分解兩層的模型和單一的SARIMA模型(除了驗證集預(yù)測的MAPE)，為最優(yōu)模型。

圖6 三種模型擬合與預(yù)測圖

表1 三種模型的擬合與預(yù)測性能指標(biāo)

討 論

本研究引入小波分析，構(gòu)建基于小波分析的SARIMA模型，利用小波對原始序列進(jìn)行分解，將可能包含噪聲的細(xì)節(jié)成分去除，只使用近似成分序列進(jìn)行建模，發(fā)現(xiàn)這種基于小波的SARIMA模型的擬合和預(yù)測性能不錯，指標(biāo)值總體上比單一的SARIMA模型降低，其性能較單一的SARIMA模型好。由于小波分解層數(shù)過多會損失可能有用的信息，所以本研究只進(jìn)行了一層和兩層分解。結(jié)果表明，只分解一層效果較好，提示分解一層后產(chǎn)生的細(xì)節(jié)成分大多包含的是隨機(jī)噪聲項。

以往手足口病時間序列預(yù)測研究多采用SARIMA模型[12-13]，其簡單實用，是一種相對比較成熟的線性時間序列預(yù)測模型。但其擬合非穩(wěn)定性時間序列的能力是有限的[14]。小波分析是一種研究非穩(wěn)定性時間序列的強(qiáng)有力的工具[15]。小波分解可以將原始序列分解為高頻成分和低頻成分，不同成分進(jìn)行不同處理[16]。通過這種方法，簡化時間序列，提高預(yù)測精度。本研究結(jié)果表明基于小波分析的SARIMA模型適合用于手足口病時間序列預(yù)測，可提高預(yù)測精度，有助于降低發(fā)病率、減輕疾病負(fù)擔(dān)。

雖然研究結(jié)果表明基于小波分析的SARIMA模型優(yōu)于單一的SARIMA模型，但手足口病影響因素比較復(fù)雜，不同地區(qū)不同時間流行模式可能不同，是否這種結(jié)合了小波的SARIMA模型能夠穩(wěn)定有效地提高手足口病預(yù)測精度尚待后續(xù)更多研究。此外，不同小波性質(zhì)不同，我們僅僅選取了Daubechies小波進(jìn)行分解，未嘗試其他小波；細(xì)節(jié)成分中或有可能包含有用的信息成分，本研究僅僅將之舍棄不用，未嘗試其他處理方法，今后研究可對其他小波和細(xì)節(jié)成分處理方式進(jìn)行更多的探索。