翟夢夢 王旭春 任 浩 全帝臣 李美晨 陳利民 仇麗霞△

【提 要】 目的 探討基于keras的LSTM模型和SARIMA模型預測我國北方省份流感樣病例數(shù)的可行性，為流感防控工作提供合理的預測方法。方法 利用國家流感中心2013-2019年北方省份的周流感監(jiān)測數(shù)據(jù)構(gòu)建LSTM模型和SARIMA模型，并進行預測。采用平均絕對誤差(MAE)、均方根誤差(RMSE)評價兩種模型的預測效果。結(jié)果 LSTM、SARIMA模型的MAE值分別為304.19、352.74，RMSE值分別為398.71、521.07；相比之下，LSTM模型的預測性能優(yōu)于SARIMA，較SARIMA模型預測性能分別提高了13.76%、23.5%。結(jié)論 基于Keras的LSTM模型的預測效果較好，優(yōu)于SARIMA模型，可為流感預測提供科學依據(jù)。

流行性感冒簡稱流感，流行期間會造成學生缺課、人員缺勤、勞動力嚴重喪失，嚴重影響人群的健康狀況[1]；流行高峰期間，發(fā)病和就診人數(shù)急劇增加，門診和住院負擔也會隨之加重[2-3]，阻礙了社會經(jīng)濟的有序發(fā)展。近年來我國流感發(fā)病數(shù)呈上升趨勢，因此，準確、合理的預測模型對政府的決策和有效防控有著重要的指導意義。

流感監(jiān)測數(shù)據(jù)具有時間序列特征，自回歸移動平均(autoregressive integrated moving average，ARIMA)模型作為傳統(tǒng)時間序列預測的代表，能夠分析數(shù)據(jù)的長期趨勢、周期性變化和隨機擾動，是疾病預測中最常用的一種方法[4]。當時間序列數(shù)據(jù)存在季節(jié)性或周期性波動時，常構(gòu)建季節(jié)性自回歸移動平均(seasonal autoregressive integrated moving average，SARIMA)模型[5]，但在實際應用中，往往需要通過差分等方法將不平穩(wěn)序列變?yōu)槠椒€(wěn)序列，這會造成信息浪費且隨著預測時間的延長，精確度隨之降低[6]。近年來，隨著計算機硬件和算法的提升，深度神經(jīng)網(wǎng)絡(deep neural network，DNN)的計算得以實現(xiàn)，與傳統(tǒng)時間序列方法相比，不需要對原始時序數(shù)據(jù)進行平穩(wěn)性處理，也能獲取數(shù)據(jù)的波動性及規(guī)律性特征[7]，其中，循環(huán)神經(jīng)網(wǎng)絡(recurrent neural network，RNN)的預測性能較好，但當時間序列較長時，RNN會出現(xiàn)梯度消失及長期記憶能力不足的問題[8-9]。1997年，Jürgen Schmidhuber等人提出了長短期記憶神經(jīng)網(wǎng)絡(long-short term memory，LSTM)，可以很好地解決上述問題，更適合處理長時間序列數(shù)據(jù)[10]。目前，LSTM已被應用于股票預測[11]、銷售預測[12]等領域，但在流感預測方面少見報道，本研究利用2013-2019年國家流感中心報告的北方省份流感樣病例建立基于keras平臺的LSTM模型，研究其對流感的預測性能，并與SARIMA模型進行比較，為我國傳染病預測提供科學依據(jù)。

資料與方法

1.資料來源

本研究使用的數(shù)據(jù)來源于國家流感中心(http：//www.chinaivdc.cn/cnic/)。國內(nèi)的流感調(diào)查研究發(fā)現(xiàn)，我國南方與北方流感流行情況存在差異，由于北方的流感樣病例報告數(shù)較南方省份有更為明顯的季節(jié)特征，更適合建立SARIMA模型，因此，整理了北方省份2013-2019年共364周的流感監(jiān)測周報數(shù)據(jù)進行分析和建模。

2.研究方法

(1)SARIMA模型

預處理過程包括：①平穩(wěn)性檢驗：根據(jù)單位根檢驗(augmented dickey-fuller test，ADF)及自相關函數(shù)圖判斷序列平穩(wěn)性。②白噪聲檢驗(Box-Ljung檢驗)：對平穩(wěn)序列進行白噪聲檢驗，若序列為白噪聲序列則分析結(jié)束，若為非白噪聲序列則進行建立SARIMA模型。

SARIMA模型表達式為：(p，d，q)(P，D，Q)S。p，q為自回歸和移動平均階數(shù)，d為差分次數(shù)，P、Q為季節(jié)自回歸和移動平均階數(shù)，D為季節(jié)差分次數(shù)。其構(gòu)建過程主要包括：①根據(jù)調(diào)整后序列的自相關圖與偏相關圖，選擇可能的p、q、Q和P值進行擬合。②運用極大似然估計法估計模型參數(shù)。③模型檢驗：最終通過最小信息量準則(the Akaike′s information criterion，AIC)和施瓦斯貝葉斯準則(Schwars′s bayesian criterion，SBC)最小原則選擇最優(yōu)模型。④利用所選的最優(yōu)模型進行預測。

(2)LSTM模型

LSTM模型于1997年被Sepp Hochreiter和 Jurgen Schmidhuber提出，后被Alex Graves、Haim Sak和Wojciech Zaremba等人逐步改進并予以應用[10]，是一種特殊的RNN網(wǎng)絡類型。LSTM模型中引入了獨特的記憶單元結(jié)構(gòu)，使其可以在學習時間序列長期依賴性的同時克服RNN的梯度消失問題，使網(wǎng)絡能夠更加有效地利用長時間序列信息，同時，LSTM模型對時間序列無平穩(wěn)性要求，可以直接對序列進行訓練，減少了人為因素的干預，使得模型具有較強的客觀性，因此，本研究選擇該模型進行分析與預測。

由于LSTM模型是RNN的一種改進模型，其模型基本結(jié)構(gòu)與RNN相同。RNN結(jié)構(gòu)如圖1所示，LSTM模型僅與RNN的隱含層模塊(單元結(jié)構(gòu)A)不同。

圖1 RNN結(jié)構(gòu)

LSTM模型是將RNN的隱含層模塊(單元結(jié)構(gòu)A)重新設置為LSTM單元結(jié)構(gòu)，其內(nèi)部結(jié)構(gòu)如圖2所示。圖2中的A分別代表一個LSTM單元。

圖2 LSTM網(wǎng)絡單元結(jié)構(gòu)

LSTM模型中的單元結(jié)構(gòu)A(LSTM單元)又稱記憶單元，由輸入門、遺忘門和輸出門組成，其中，a表示遺忘門，決定傳來的信息是否要丟棄；b表示輸入門，與c共同更新當前時刻的單元狀態(tài)，具體是通過將it和mt相乘更新當前時刻的單元狀態(tài)Ct，并將其傳遞到下一個記憶單元；d表示輸出門，決定當前時刻輸出的值，由ot和經(jīng)tanh函數(shù)處理的Ct共同決定。模型通過這三種門控結(jié)構(gòu)控制信息在記憶單元和神經(jīng)網(wǎng)絡中的流動[13]。ft、it、ot分別代表三個門結(jié)構(gòu)狀態(tài)，mt代表單元輸入狀態(tài)，ht-1代表t-1時刻的輸出、xt代表t時刻的輸入。具體計算公式如下[14]：

ft=σ(Wfhht-1+Wfxxt+bf)

it=σ(Wihht-1+Wixxt+bi)

mt=tanh(Wmhht-1+Wmxxt+bm)

Ct=ftCt-1+immt

ot=σ(Wohht-1+Woxxt+bo)

ht=ottanh(Ct)

式中，W表示連接兩層的權(quán)重矩陣，σ表示sigmoid激活函數(shù)，b是對應的偏置項。

LSTM模型構(gòu)建過程：①將原始數(shù)據(jù)進行預處理，包括數(shù)據(jù)歸一化及數(shù)據(jù)重構(gòu)；②轉(zhuǎn)換為3D數(shù)據(jù)：將預處理過的數(shù)據(jù)轉(zhuǎn)換為適合LSTM要求的3D數(shù)據(jù)；③基于Python的深度學習框架Tensorflow下的keras平臺來建立LSTM模型，使用訓練集數(shù)據(jù)對模型進行訓練；④使用訓練后的模型在測試集上進行；⑤預測數(shù)據(jù)反歸一化。

數(shù)據(jù)歸一化：為了提高模型的訓練速度和預測精度[15]，本文使用min-max normalization方法對原始數(shù)據(jù)進行歸一化處理。

式中：X*為數(shù)據(jù)歸一化后的值；X為原始數(shù)據(jù)，Xmax和Xmin分別為其最大值和最小值。

數(shù)據(jù)重構(gòu)：由于流感數(shù)據(jù)為一維時間序列，而LSTM模型的數(shù)據(jù)樣本要有輸入X和對應的輸出Y，因此，本文通過滑窗法進行數(shù)據(jù)重構(gòu)。

滑窗法是使用過去時間點觀測值來預測下一時間點預測值的方法，具體重構(gòu)過程如圖3所示，即假設原始時間序列為具有m個時間點的觀測值序列，將t時刻之前的k個時間點的觀測值作為輸入Xt，t時刻的觀測值作為輸出Yt。以此類推，將原始時間序列重構(gòu)為有輸入和輸出的數(shù)據(jù)集。圖3中將k設置為6(即使用t時刻前6個時間點的觀測值作為輸入)，此時Xt為[[3]，[4]，[5]，[6]，[7]，[8]]，是一個2D張量，t時刻的Yt為[[9]]，同樣是2D張量。Xt+1為[[4]，[5]，[6]，[7]，[8]，[9]]，對應的Yt+1為[[10]]。以此類推，則重構(gòu)后的數(shù)據(jù)集包含m-k個樣本點。其中，1，2，3，4…8，9，10…m分別代表各時間點觀測值歸一化后的數(shù)值。

圖3 滑窗法數(shù)據(jù)重構(gòu)示意圖

數(shù)據(jù)轉(zhuǎn)換：由于LSTM模型要求輸入數(shù)據(jù)的格式為(samples，time_steps，input_dim)的3D張量(3D數(shù)據(jù))。其中，samples代表樣本量，time_steps代表時間步長，input_dim代表每個時間步的數(shù)據(jù)維度。經(jīng)數(shù)據(jù)重構(gòu)后(圖3)，得到的輸入樣本Xt-2、Xt-1、Xt、Xt+1、Xt+2…Xm-6均為2D張量，將重構(gòu)后的輸入樣本按順序組合在一起時就轉(zhuǎn)換成了LSTM所要求的3D張量。由于時間步長不同，重構(gòu)后的樣本量也不相同，當時間步長為k時，樣本量為m-k，且流感序列為一維數(shù)據(jù)，因此，LSTM模型的輸入數(shù)據(jù)格式為(m-k，k，1)。

(3)模型效果評價

采用平均絕對誤差(mean absolute error，MAE)、均方根誤差(root mean squared absolute error，RMSE)指標來評價模型的預測性能。

(1)

(2)

3.統(tǒng)計分析

采用SAS 9.2軟件建立SARIMA模型，采用python 3.7.1進行單位根檢驗、白噪聲檢驗、自相關圖和偏自相關圖；采用Pycharm的keras軟件包實現(xiàn)LSTM模型。

結(jié) 果

1.流感樣病例時間分布

2013-2019年我國北方省份流感流行概況如圖4、5所示?？梢娏鞲袠硬±磕瓿试鲩L趨勢，且每年的冬春交替時期是流感發(fā)病的高峰期。

圖4 2013-2019年北方省份流感樣病例時序圖

2.SARIMA模型

由圖4、圖5可知該序列呈每年遞增的趨勢和明顯的季節(jié)性特征，原始序列為不平穩(wěn)序列(圖6)。將2013-2019年共364周數(shù)據(jù)的前4/5設置為訓練集，后1/5設置為測試集。進行1階差分及52步差分，差分后序列平穩(wěn)且不能視為白噪聲(表1)，季節(jié)性周期為52周(即S=52)，初步選定SARIMA(p，1，q)(P，1，Q)52模型。

圖5 北方省份年度流感樣病例數(shù)

圖6 原始序列自相關與偏相關圖

表1 差分后序列平穩(wěn)性及白噪聲檢驗

根據(jù)文獻顯示，參數(shù)超過2階的情況極為少見，且通過差分后序列的自相關和偏相關圖(圖7)確定可能的p、q、P和Q值，確定5個備選模型(表2)，根據(jù)AIC，SBC最小原則選擇最優(yōu)模型，最終確定模型為SARIMA(0，1，1)(0，1，1)52。該模型的殘差白噪聲檢驗結(jié)果為χ2=17.124，P=0.378，表明其殘差序列為白噪聲。

表2 備選模型參數(shù)估計及參數(shù)檢驗

圖7 差分后序列自相關與偏相關圖

從QQ圖看出殘差服從正態(tài)分布，殘差序列為白噪聲，再次表明所構(gòu)建的模型是有效的。利用SARIMA(0，1，1)(0，1，1)52模型預測我國北方省份2019年1月-12月的流感樣病例數(shù)，結(jié)果如圖9所示。

圖9 SARIMA(0，1，1)(0，1，1)52模型預測圖

3.基于keras平臺的LSTM模型

由于LSTM模型對時間序列無平穩(wěn)性要求，可以直接建立模型。同樣地將前4/5(291)的流感數(shù)據(jù)設置為訓練集，后1/5(73)設置為測試集，建立LSTM模型。

圖8 SARIMA模型殘差QQ圖

本文的隱藏層設置為單隱層，節(jié)點數(shù)設置為4，輸出層設置為1，采用Adam算法優(yōu)化參數(shù)，其他為默認參數(shù)。由于時間步長會影響待預測樣本的長度，原始數(shù)據(jù)經(jīng)歸一化及數(shù)據(jù)重構(gòu)后將其轉(zhuǎn)換為LSTM模型輸入要求的3D數(shù)據(jù)格式，經(jīng)數(shù)據(jù)轉(zhuǎn)化可知：其數(shù)據(jù)格式為(m-k，k，1)，其中，m為原始時間序列的樣本數(shù)，k為時間步長。經(jīng)過數(shù)據(jù)重構(gòu)，時間步為1的LSTM模型訓練集和測試集樣本數(shù)為290和72，時間步為52的LSTM模型樣本數(shù)為239和21。不同時間步長下的輸入數(shù)據(jù)格式如表3所示。

表3 不同時間步長下的輸入數(shù)據(jù)格式

當模型的時間步設置為1時(即將前1周數(shù)據(jù)作為輸入，后1周數(shù)據(jù)作為輸出)，結(jié)果見圖10。從圖中可以看出，在數(shù)據(jù)不做差分等處理的前提下，LSTM模型抓住了數(shù)據(jù)線性增長和波動率逐漸增加的兩大趨勢，這是經(jīng)典的時間序列分析模型不容易做到的；但由于原始數(shù)據(jù)存在周期性的結(jié)構(gòu)特點，將時間步設置為52(即將前52周數(shù)據(jù)作為輸入，后1周數(shù)據(jù)作為輸出)進行預測，結(jié)果見圖11。并與時間步為1的模型進行比較。

圖10 LSTM(1)模型預測圖及100次預測丟失率圖

圖11 LSTM(52)模型預測圖及100次預測丟失率圖

時間步為1的LSTM模型預測丟失率穩(wěn)定在0.0096，高于時間步為52的LSTM模型預測丟失率穩(wěn)定在0.0039，其中，丟失率(loss)代表訓練集樣本的預測值與真實值的差異程度，用于評價LSTM模型的訓練效果。通過比較兩模型的丟失率，表明所構(gòu)建的時間步為52的LSTM模型較好。

4.模型對比分析

運用SARIMA模型和LSTM模型對北方省份2019年1-12月流感樣病例數(shù)進行預測，結(jié)果顯示：SARIMA模型、LSTM(1)、LSTM(52)模型對于流感的預測趨勢與真實流感樣病例數(shù)趨勢基本一致，說明SARIMA、LSTM模型都能較好預測流感發(fā)病趨勢。SARIMA模型預測效果居中，其MAE、RMSE分別為352.74、521.07；LSTM(1)預測效果最差，其MAE、RMSE分別為448.79、784.09，較SARIMA模型分別增加96.05、263.02，預測性能分別降低了27.23%、50.50%；LSTM(52)預測效果最好，其MAE、RMSE分別為304.19、398.71，較SARIMA模型分別減少48.55、122.36，預測性能分別提高了13.76%、23.50%；表明對于具有周期性特征的時間序列而言，考慮周期特征的LSTM模型在流感預測性能最好(表4)。

表4 三種模型預測性能對比

討 論

由于流感病毒極易發(fā)生變異，造成大流行[16]，加重衛(wèi)生服務負擔。因此，及時了解流感流行趨勢、早期發(fā)現(xiàn)疫情是預防控制的關鍵。SARIMA模型著重分析流感樣病例數(shù)隨時間的變化規(guī)律，可以將其他未知的影響因素綜合于時間序列中，最終通過分析流感樣病例數(shù)之間的自相關性和依存關系來進行預測，SARIMA模型要求序列為平穩(wěn)的非白噪聲序列，對于非平穩(wěn)的時間序列，需要對序列進行差分達到平穩(wěn)后進行分析，而LSTM模型則可以直接對時間序列進行訓練并預測，相比于傳統(tǒng)的統(tǒng)計學分析模型更具優(yōu)勢[17]，因此本文建立SARIMA模型和LSTM模型，并根據(jù)時間周期性特征對LSTM模型進行調(diào)參，采用MAE和RMSE兩個指標來比較兩模型在預測流感樣病例數(shù)的準確性。

結(jié)果顯示，3個模型對我國北方省份流感的預測趨勢與真實流感病例數(shù)趨勢基本一致，但相比較來說，LSTM(52)，即考慮時間周期性，模型的預測效果最好，優(yōu)于SARIMA模型和LSTM(1)模型?？赡苁怯捎贚STM(1)模型只用前1周的流感發(fā)病人數(shù)預測下一周，沒有考慮流感發(fā)病的周期性，而LSTM(52)模型克服了平移錯位的現(xiàn)象，提高了預測精度，這也與李琳等[8]利用LSTM模型分析新疆地區(qū)慢性阻塞性肺病的月門診量變化趨勢研究結(jié)果相一致。因此，與傳統(tǒng)的時間序列預測模型SARIMA相比，本文建立的基于keras平臺的LSTM模型可以對未來一段時間我國流感樣病例數(shù)進行預測，為流感有效防控及相關部門進行資源優(yōu)化配置提供一定的參考。