SARIMA模型和LSTM神經(jīng)網(wǎng)絡(luò)在預(yù)測深圳市寶安區(qū)手足口病疫情中的應(yīng)用

陳春艷，陳億雄，趙執(zhí)揚(yáng)，李 靜，李淑珍*

(1山西醫(yī)科大學(xué)公共衛(wèi)生學(xué)院流行病學(xué)教研室，太原 030001；2深圳市寶安區(qū)疾病預(yù)防控制中心傳防科；*通訊作者,E-mail:lszlym@163.com)

手足口病(hand-foot-mouth disease, HFMD)是由多種腸道病毒感染引起的常見急性傳染病，常見于5歲以下兒童[1]。臨床上以手、足和口腔等部位出現(xiàn)皰疹、斑丘疹為主要癥狀[2]。多數(shù)患者癥狀較輕，可自愈，少數(shù)重癥患者可出現(xiàn)中樞神經(jīng)系統(tǒng)和/或心血管系統(tǒng)并發(fā)癥如肺循環(huán)衰竭、腦膜炎等[3]，甚至死亡。手足口病傳播途徑復(fù)雜、防控難度大，我國于2008年將其列為法定報(bào)告的丙類傳染病[4]。

時間序列分析[5]是將某種現(xiàn)象某一個統(tǒng)計(jì)指標(biāo)在不同時間上的各個數(shù)值，按時間先后順序排列而形成的序列，旨在預(yù)測未來事件發(fā)展的趨勢和規(guī)律，現(xiàn)已被各領(lǐng)域廣泛應(yīng)用。時間序列分析中的季節(jié)性差分自回歸移動平均模型(seasonal autoregressive integrated moving average model, SARIMA)能夠分析發(fā)病數(shù)據(jù)的趨勢性、季節(jié)性、周期性以及隨機(jī)性的波動，最常用于預(yù)測疾病的發(fā)生及發(fā)展規(guī)律[6,7]，但SARIMA模型要求其時間序列服從平穩(wěn)性，當(dāng)序列不平穩(wěn)時往往需要通過差分等方法將其變?yōu)槠椒€(wěn)序列，因此會損失一定的信息造成預(yù)測準(zhǔn)確性降低。近年來，隨著深度學(xué)習(xí)理論的出現(xiàn)和數(shù)值計(jì)算能力的提升，基于循環(huán)神經(jīng)網(wǎng)絡(luò)(recurrent neural networks, RNN)改進(jìn)的一種算法即長短時記憶神經(jīng)網(wǎng)絡(luò)(long short term memory, LSTM)逐漸成為時間序列的分析方法之一[8,9]。

廣東省深圳市是手足口病的高發(fā)城市[10,11]。寶安區(qū)是深圳市人口數(shù)最多的區(qū)，達(dá)447萬人[12]，是深圳市手足口病的高發(fā)地區(qū)[13]。因?yàn)槭苄鹿谝咔榈挠绊懀?020—2021年通過“中國疾病預(yù)防控制信息系統(tǒng)”報(bào)告的深圳市寶安區(qū)的HFMD發(fā)病數(shù)與往年相比，數(shù)據(jù)質(zhì)量穩(wěn)定性受到影響，容易導(dǎo)致構(gòu)建的模型具有較大的預(yù)測誤差。為減少誤差，本文選取2009年1月至2019年12月深圳市寶安區(qū)的HFMD發(fā)病數(shù)據(jù)資料進(jìn)行研究。通過應(yīng)用SARIMA模型和LSTM神經(jīng)網(wǎng)絡(luò)構(gòu)建深圳市寶安區(qū)手足口病時間序列模型，預(yù)測發(fā)病趨勢，為今后HFMD的預(yù)防和控制提供理論依據(jù)。

1 資料與方法

1.1 一般資料

2009年1月至2019年12月深圳市寶安區(qū)人口數(shù)資料來源于深圳市寶安區(qū)統(tǒng)計(jì)局，2009年1月至2019年12月深圳市寶安區(qū)手足口病發(fā)病數(shù)據(jù)來源于“中國疾病預(yù)防控制信息系統(tǒng)”。本研究以2009年1月至2018年12月的HFMD月發(fā)病率作為訓(xùn)練集分別構(gòu)建SARIMA模型和LSTM神經(jīng)網(wǎng)絡(luò)，預(yù)測2019年1—12月的HFMD月發(fā)病率。

1.2 研究方法

1.2.1 SARIMA模型 由于手足口病具有季節(jié)性波動，本研究構(gòu)建季節(jié)性差分自回歸移動平均模型(seasonal autoregressive integrated moving average model, SARIMA)，即SARIMA(p，d，q)(P，D，Q)s。構(gòu)建SARIMA模型的步驟為：首先進(jìn)行自相關(guān)圖檢驗(yàn)、單位根檢驗(yàn)法(augmented dickey-fullert，ADF)判斷數(shù)據(jù)是否平穩(wěn)(P<0.05代表序列平穩(wěn))，若序列不平穩(wěn)則需差分直至成為平穩(wěn)序列。隨后對該序列進(jìn)行白噪聲檢驗(yàn)即Ljung-Box檢驗(yàn)，當(dāng)序列成為非白噪聲時，則可構(gòu)建SARIMA模型。本研究采用Python網(wǎng)格搜索來自動擬合SARIMA模型，根據(jù)赤池信息(AIC最小)準(zhǔn)則，采用條件最小二乘法估計(jì)模型參數(shù)，并對模型參數(shù)進(jìn)行統(tǒng)計(jì)學(xué)檢驗(yàn)，選擇相對最優(yōu)模型，通過模型殘差白噪聲檢驗(yàn)判斷模型是否擬合成功(Ljung-Box檢驗(yàn)中P>0.05，表示殘差為白噪聲)。

1.2.2 LSTM神經(jīng)網(wǎng)絡(luò) LSTM神經(jīng)網(wǎng)絡(luò)是由Hochreiter提出并由Graves改進(jìn)的一種常見的循環(huán)神經(jīng)網(wǎng)絡(luò)[14]。一個LSTM單元結(jié)構(gòu)的核心在于它的神經(jīng)元狀態(tài)(cell state)。LSTM包括3個門控結(jié)構(gòu)[8]，即“遺忘門”“輸入門”和“輸出門”。這3種門控結(jié)構(gòu)可選擇性地控制信息通過，信息傳遞順序?yàn)?先輸入ht-1、xt和Ct-1，根據(jù)sigmoid、tanh函數(shù)和下列公式[8]計(jì)算ft、it、ot，其中ht-1代表t-1時刻的輸出值，xt代表t時刻的輸入值，Ct-1代表t-1時刻的單元狀態(tài)，ft、it、ot分別代表遺忘門狀態(tài)、輸入門狀態(tài)、輸出門狀態(tài)[15]。

遺忘門狀態(tài)：ft=σ(Wf·[Ct-1,ht-1,xt]+bf)

輸入門狀態(tài)：it=σ(Wi·[Ct-1,ht-1,xt]+bi)

輸出門狀態(tài)：ot=σ(Wo·[Ct-1,ht-1,xt]+bo)

單元輸入狀態(tài):

隱藏層輸出狀態(tài)：ht=ot×tanh(Ct)

其中W是權(quán)重矩陣，b是偏倚項(xiàng)，σ表示sigmoid函數(shù)。

1.3 模型性能評價(jià)

運(yùn)用均方誤差(mean squared error, MSE)、均方根誤差(root mean squared error, RMSE)、平均絕對誤差(mean absolute error, MAE)、平均絕對百分比誤差(mean absolute percentage error, MAPE)比較兩種模型對深圳市寶安區(qū)HFMD發(fā)病趨勢的擬合及預(yù)測效果。

1.4 統(tǒng)計(jì)分析

采用SPSS 22.0軟件對數(shù)據(jù)進(jìn)行描述性分析，采用χ2檢驗(yàn)進(jìn)行率的比較，檢驗(yàn)水準(zhǔn)ɑ=0.05；采用python3.10軟件分別構(gòu)建SARIMA模型和LSTM神經(jīng)網(wǎng)絡(luò)。

2 結(jié)果

2.1 發(fā)病率分析

2.1.1 年發(fā)病率分析 2009—2019年深圳市寶安區(qū)累計(jì)報(bào)告手足口病82 632例，年均發(fā)病率為291.93/10萬，2017年(575.34/10萬)及2018年(699.84/10萬)的HFMD發(fā)病率較高(見表1)。

2.1.2 月發(fā)病率分析 除2009—2012年外，其他年份的HFMD月發(fā)病率均呈現(xiàn)典型的雙峰模式。病例從每年的3月開始增多，到第4—7月達(dá)到第1個峰值，隨后第9—11月達(dá)到第2個峰值(見圖1)，即多發(fā)于夏秋季節(jié)，夏季較秋季高發(fā)，具有明顯的季節(jié)性。

表1 2009—2019年深圳市寶安區(qū)手足口病發(fā)病的分布Table 1 Distribution of incidence of HFMD in Bao’an district of Shenzhen city in 2009—2019

圖1 2009—2019年深圳市寶安區(qū)手足口病月發(fā)病率分布Figure 1 Distribution of monthly incidence of HFMD in Bao’an district of Shenzhen city in 2009—2019

2.2 SARIMA模型

2.2.1 平穩(wěn)性檢驗(yàn)和白噪聲檢驗(yàn) 對原始序列依次進(jìn)行平穩(wěn)性檢驗(yàn)和白噪聲檢驗(yàn)，經(jīng)ADF檢驗(yàn)發(fā)現(xiàn)，序列存在單位根(χ2=-0.018，P>0.05，見表2)，對原始序列的自相關(guān)圖、偏自相關(guān)圖進(jìn)行分析，自相關(guān)系數(shù)超過2倍標(biāo)準(zhǔn)差(見圖2)，綜合考慮上述結(jié)果，認(rèn)為原始序列為不平穩(wěn)序列。對序列進(jìn)行季節(jié)性差分得到平穩(wěn)序列(χ2=-3.091，P<0.05)，隨后對該序列進(jìn)行Ljung-Box檢驗(yàn)，序列滿足非白噪聲要求(χ2=15.841，P<0.05，見表2)。

2.2.2 模型優(yōu)化 根據(jù)序列的差分次數(shù)確定非季節(jié)性與季節(jié)性差分的值，即d為0，D為1，初步選用SARIMA(p，0，q)(P，1，Q)12模型。根據(jù)既往文獻(xiàn)經(jīng)驗(yàn)[16]p、q、P及Q取值范圍為0～2，采用Python網(wǎng)格搜索自動擬合SARIMA模型，根據(jù)AIC最小原則和解釋變量盡可能有意義，最終選擇SARIMA(0，0，2)(0，1，2)12為相對最優(yōu)模型(AIC=752.094，BIC=764.066)，對模型的殘差序列進(jìn)行Ljung-Box檢驗(yàn)，結(jié)果為白噪聲(P=0.510，見表3)，說明該模型數(shù)據(jù)提取完全，效果滿意。

表2 序列的平穩(wěn)性及白噪聲檢驗(yàn)Table 2 Stationarity of sequences and Ljung-Box test

圖2 原始序列及差分序列的自相關(guān)圖與偏自相關(guān)圖Figure 2 ACF and PACF of original sequence and difference sequence

表3 SARIMA(0，0，2)(0，1，2)12模型參數(shù)估計(jì)和擬合優(yōu)度統(tǒng)計(jì)量結(jié)果Table 3 Parameter estimation and goodness-of-fit statistics for SARIMA(0，0，2)(0，1，2)12 models

2.3 LSTM神經(jīng)網(wǎng)絡(luò)

因?yàn)槭肿憧诓〉臄?shù)據(jù)周期為12個月，因此LSTM神經(jīng)網(wǎng)絡(luò)的窗口長度設(shè)置為12，即輸入層節(jié)點(diǎn)數(shù)為12，預(yù)測下個月手足口病發(fā)病率，輸出層節(jié)點(diǎn)數(shù)為1。

本文在單隱層的結(jié)構(gòu)下以隱藏層節(jié)點(diǎn)數(shù)為2的冪次方進(jìn)行試驗(yàn)，當(dāng)隱藏層節(jié)點(diǎn)數(shù)為1 024時，模型的RMSE、MAE、MSE 3個評價(jià)指標(biāo)均最小(見表4)。

固定隱藏層節(jié)點(diǎn)數(shù)為1 024，設(shè)置隱藏層層數(shù)為1～5逐個進(jìn)行試驗(yàn)，結(jié)果見表5。當(dāng)隱藏層層數(shù)設(shè)置為1時，LSTM神經(jīng)網(wǎng)絡(luò)的誤差值均最低。

綜上，本文設(shè)置時間步長為12、隱藏層層數(shù)為1、隱藏層節(jié)點(diǎn)數(shù)為1 024、迭代次數(shù)為50時，構(gòu)建LSTM神經(jīng)網(wǎng)絡(luò)。

表4 隱藏層節(jié)點(diǎn)數(shù)對模型預(yù)測性能的影響Table 4 The influence of the number of hidden layer nodes on the prediction performance of the model

2.4 兩種模型擬合及預(yù)測深圳市寶安區(qū)HFMD發(fā)病的效果評價(jià)

在本研究中，采用SARIMA(0，0，2)(0，1，2)12和LSTM神經(jīng)網(wǎng)絡(luò)分別對2010年1月至2018年12月深圳市寶安區(qū)的HFMD月發(fā)病率進(jìn)行擬合。兩種模型對2019年1—12月深圳市寶安區(qū)的HFMD月發(fā)病率進(jìn)行預(yù)測，SARIMA的預(yù)測值出現(xiàn)了雙峰分布，而LSTM的預(yù)測值呈單峰分布且與實(shí)際月發(fā)病率基本吻合(見圖3)。

表5 隱藏層層數(shù)對模型預(yù)測性能的影響Table 5 The influence of the number of hidden layers on the prediction performance of the model

為客觀評價(jià)模型擬合及預(yù)測性能，使用MSE、RMSE、MAE及MAPE對兩模型進(jìn)行比較，在擬合性能中LSTM模型的MSE、RMSE均高于SARIMA模型，而MAE、MAPE均低于SARIMA模型(見表6)，表明兩種模型的擬合性能基本一致；而在預(yù)測性能中LSTM模型的MSE、RMSE、MAE及MAPE 4個指標(biāo)均低于SARIMA模型。因此，LSTM神經(jīng)網(wǎng)絡(luò)的預(yù)測性能高于SARIMA模型，表明LSTM神經(jīng)網(wǎng)絡(luò)能更好地預(yù)測深圳市寶安區(qū)HFMD發(fā)病趨勢。

黑色虛線左側(cè)為兩種模型的擬合效果，右側(cè)為兩種模型的預(yù)測效果圖3 SARIMA和LSTM模型對2010年1月至2018年12月深圳市寶安區(qū)的HFMD月發(fā)病率擬合及預(yù)測效果的比較Figure 3 Comparison of the performance of fitting and predicting monthly incidence rate of HFMD in Bao’an district of Shenzhen city in 2010—2018 between SARIMA model and LSTM model

表6 兩種模型的擬合及預(yù)測性能的誤差值對比Table 6 Comparison of fitting and prediction errors between the two models

3 討論

自2008年我國將手足口病納入丙類傳染病管理以來，手足口病的發(fā)病率一直位于我國丙類傳染病的首位[17]。手足口病是由不同腸道病毒感染所致，不同病原體之間存在差異，易在托幼機(jī)構(gòu)及學(xué)校呈暴發(fā)流行，這類聚集性疫情的發(fā)生給國家和地區(qū)的衛(wèi)生保健系統(tǒng)造成較大的經(jīng)濟(jì)負(fù)擔(dān)[18]。本研究中，深圳市寶安區(qū)2014—2019年手足口病發(fā)病率顯著高于2009—2013年，可能原因跟該地區(qū)傳染病疫情監(jiān)測系統(tǒng)逐漸完善、外來人口數(shù)增多、人口流動性強(qiáng)等密切相關(guān)。因此，了解手足口病流行趨勢，構(gòu)建恰當(dāng)?shù)氖肿憧诓☆A(yù)測模型，可為相關(guān)部門制定手足口病防控策略提供科學(xué)依據(jù)。

本研究構(gòu)建SARIMA模型時，先對原始序列進(jìn)行平穩(wěn)性及白噪聲檢驗(yàn)，判斷原始序列是否滿足平穩(wěn)非白噪聲，其次觀察序列的ACF、PACF圖確定模型的參數(shù)范圍，通過Python網(wǎng)格搜索自動擬合較佳參數(shù)，最后運(yùn)用AIC最小原則，選取了SARIMA(0，0，2)(0，1，2)12為相對最優(yōu)模型。用該模型分別對深圳市寶安區(qū)的HFMD月發(fā)病率進(jìn)行擬合及預(yù)測時，誤差指數(shù)較小，擬合及預(yù)測原始序列的趨勢性和周期性的效果較好，表明SARIMA模型預(yù)測深圳市寶安區(qū)HFMD發(fā)病趨勢的效果尚佳，與韓玲等[19]使用河北省HFMD發(fā)病數(shù)建立SARIMA(1，1，1)(0，1，1)12模型、劉濤等[20]使用山東省HFMD發(fā)病數(shù)建立SARIMA(1，0，1)(0，1，0)12模型的研究結(jié)果一致。

本研究對LSTM神經(jīng)網(wǎng)絡(luò)的隱藏層層數(shù)及節(jié)點(diǎn)數(shù)進(jìn)行了參數(shù)的調(diào)整，設(shè)置時間步長為12、隱藏層層數(shù)為1、隱藏層節(jié)點(diǎn)數(shù)為1 024、迭代次數(shù)為50時構(gòu)建的LSTM神經(jīng)網(wǎng)絡(luò)為較優(yōu)模型，該模型對深圳市寶安區(qū)HFMD月發(fā)病率的擬合及預(yù)測性能較好，與馮一平等[21]研究認(rèn)為LSTM神經(jīng)網(wǎng)絡(luò)能較好地預(yù)測HFMD發(fā)病一致。但是，選擇最佳參數(shù)構(gòu)建最優(yōu)LSTM神經(jīng)網(wǎng)絡(luò)較難實(shí)現(xiàn)，即使部分參數(shù)可以根據(jù)數(shù)據(jù)及疾病特征進(jìn)行大致估計(jì)，但構(gòu)建最佳模型的幾率仍然較小，后續(xù)將進(jìn)一步深入研究該問題。

本文中SARIMA模型預(yù)測性能的MSE、RMSE、MAE、MAPE均大于LSTM神經(jīng)網(wǎng)絡(luò)，可能原因[22]在于構(gòu)建SARIMA模型的前提是要求原始序列平穩(wěn)，當(dāng)原始序列不平穩(wěn)時，需要進(jìn)行差分直至成為平穩(wěn)非白噪聲序列，導(dǎo)致SARIMA模型對HFMD發(fā)病數(shù)據(jù)提取不充分，故造成信息利用不充分而產(chǎn)生預(yù)測誤差，而LSTM神經(jīng)網(wǎng)絡(luò)對數(shù)據(jù)本身是否平穩(wěn)并無要求，并且其對長期的序列具有較好的預(yù)測效果[23]。本研究表明LSTM神經(jīng)網(wǎng)絡(luò)能更好地預(yù)測深圳市寶安區(qū)HFMD發(fā)病趨勢，與高秋菊等[24]預(yù)測石家莊市HFMD發(fā)病趨勢的研究結(jié)果一致。但本文存在一定的局限性，手足口病發(fā)病易受多種因素影響，本文僅使用了手足口病監(jiān)測資料，在后續(xù)構(gòu)建模型過程中應(yīng)納入影響手足口病的其他因素，進(jìn)一步提高模型的準(zhǔn)確度，進(jìn)而提高手足口病發(fā)病預(yù)測的效果，為手足口病的防控提供更可靠的理論依據(jù)。