SARIMA模型和LSTM神經網絡在預測深圳市寶安區手足口病疫情中的應用

陳春艷，陳億雄，趙執揚，李 靜，李淑珍*

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

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

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

廣東省深圳市是手足口病的高發城市[10,11]。寶安區是深圳市人口數最多的區，達447萬人[12]，是深圳市手足口病的高發地區[13]。因為受新冠疫情的影響，2020—2021年通過“中國疾病預防控制信息系統”報告的深圳市寶安區的HFMD發病數與往年相比，數據質量穩定性受到影響，容易導致構建的模型具有較大的預測誤差。為減少誤差，本文選取2009年1月至2019年12月深圳市寶安區的HFMD發病數據資料進行研究。通過應用SARIMA模型和LSTM神經網絡構建深圳市寶安區手足口病時間序列模型，預測發病趨勢，為今后HFMD的預防和控制提供理論依據。

1 資料與方法

1.1 一般資料

2009年1月至2019年12月深圳市寶安區人口數資料來源于深圳市寶安區統計局，2009年1月至2019年12月深圳市寶安區手足口病發病數據來源于“中國疾病預防控制信息系統”。本研究以2009年1月至2018年12月的HFMD月發病率作為訓練集分別構建SARIMA模型和LSTM神經網絡，預測2019年1—12月的HFMD月發病率。

1.2 研究方法

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

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

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

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

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

單元輸入狀態:

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

其中W是權重矩陣，b是偏倚項，σ表示sigmoid函數。

1.3 模型性能評價

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

1.4 統計分析

采用SPSS 22.0軟件對數據進行描述性分析，采用χ2檢驗進行率的比較，檢驗水準ɑ=0.05；采用python3.10軟件分別構建SARIMA模型和LSTM神經網絡。

2 結果

2.1 發病率分析

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

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

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

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

2.2 SARIMA模型

2.2.1 平穩性檢驗和白噪聲檢驗 對原始序列依次進行平穩性檢驗和白噪聲檢驗，經ADF檢驗發現，序列存在單位根(χ2=-0.018，P>0.05，見表2)，對原始序列的自相關圖、偏自相關圖進行分析，自相關系數超過2倍標準差(見圖2)，綜合考慮上述結果，認為原始序列為不平穩序列。對序列進行季節性差分得到平穩序列(χ2=-3.091，P<0.05)，隨后對該序列進行Ljung-Box檢驗，序列滿足非白噪聲要求(χ2=15.841，P<0.05，見表2)。

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

表2 序列的平穩性及白噪聲檢驗Table 2 Stationarity of sequences and Ljung-Box test

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

表3 SARIMA(0，0，2)(0，1，2)12模型參數估計和擬合優度統計量結果Table 3 Parameter estimation and goodness-of-fit statistics for SARIMA(0，0，2)(0，1，2)12 models

2.3 LSTM神經網絡

因為手足口病的數據周期為12個月，因此LSTM神經網絡的窗口長度設置為12，即輸入層節點數為12，預測下個月手足口病發病率，輸出層節點數為1。

本文在單隱層的結構下以隱藏層節點數為2的冪次方進行試驗，當隱藏層節點數為1 024時，模型的RMSE、MAE、MSE 3個評價指標均最小(見表4)。

固定隱藏層節點數為1 024，設置隱藏層層數為1～5逐個進行試驗，結果見表5。當隱藏層層數設置為1時，LSTM神經網絡的誤差值均最低。

綜上，本文設置時間步長為12、隱藏層層數為1、隱藏層節點數為1 024、迭代次數為50時，構建LSTM神經網絡。

表4 隱藏層節點數對模型預測性能的影響Table 4 The influence of the number of hidden layer nodes on the prediction performance of the model

2.4 兩種模型擬合及預測深圳市寶安區HFMD發病的效果評價

在本研究中，采用SARIMA(0，0，2)(0，1，2)12和LSTM神經網絡分別對2010年1月至2018年12月深圳市寶安區的HFMD月發病率進行擬合。兩種模型對2019年1—12月深圳市寶安區的HFMD月發病率進行預測，SARIMA的預測值出現了雙峰分布，而LSTM的預測值呈單峰分布且與實際月發病率基本吻合(見圖3)。

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

為客觀評價模型擬合及預測性能，使用MSE、RMSE、MAE及MAPE對兩模型進行比較，在擬合性能中LSTM模型的MSE、RMSE均高于SARIMA模型，而MAE、MAPE均低于SARIMA模型(見表6)，表明兩種模型的擬合性能基本一致；而在預測性能中LSTM模型的MSE、RMSE、MAE及MAPE 4個指標均低于SARIMA模型。因此，LSTM神經網絡的預測性能高于SARIMA模型，表明LSTM神經網絡能更好地預測深圳市寶安區HFMD發病趨勢。

黑色虛線左側為兩種模型的擬合效果，右側為兩種模型的預測效果圖3 SARIMA和LSTM模型對2010年1月至2018年12月深圳市寶安區的HFMD月發病率擬合及預測效果的比較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 兩種模型的擬合及預測性能的誤差值對比Table 6 Comparison of fitting and prediction errors between the two models

3 討論

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

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

本研究對LSTM神經網絡的隱藏層層數及節點數進行了參數的調整，設置時間步長為12、隱藏層層數為1、隱藏層節點數為1 024、迭代次數為50時構建的LSTM神經網絡為較優模型，該模型對深圳市寶安區HFMD月發病率的擬合及預測性能較好，與馮一平等[21]研究認為LSTM神經網絡能較好地預測HFMD發病一致。但是，選擇最佳參數構建最優LSTM神經網絡較難實現，即使部分參數可以根據數據及疾病特征進行大致估計，但構建最佳模型的幾率仍然較小，后續將進一步深入研究該問題。

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