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

(1山西醫(yī)科大學公共衛(wèi)生學院流行病學教研室，太原 030001；2深圳市寶安區(qū)疾病預防控制中心傳防科；*通訊作者,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]，甚至死亡。手足口病傳播途徑復雜、防控難度大，我國于2008年將其列為法定報告的丙類傳染病[4]。

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

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

1 資料與方法

1.1 一般資料

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

1.2 研究方法

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

1.2.2 LSTM神經(jīng)網(wǎng)絡 LSTM神經(jīng)網(wǎng)絡是由Hochreiter提出并由Graves改進的一種常見的循環(huán)神經(jīng)網(wǎng)絡[14]。一個LSTM單元結構的核心在于它的神經(jīng)元狀態(tài)(cell state)。LSTM包括3個門控結構[8]，即“遺忘門”“輸入門”和“輸出門”。這3種門控結構可選擇性地控制信息通過，信息傳遞順序為:先輸入ht-1、xt和Ct-1，根據(jù)sigmoid、tanh函數(shù)和下列公式[8]計算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是權重矩陣，b是偏倚項，σ表示sigmoid函數(shù)。

1.3 模型性能評價

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

1.4 統(tǒng)計分析

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

2 結果

2.1 發(fā)病率分析

2.1.1 年發(fā)病率分析 2009—2019年深圳市寶安區(qū)累計報告手足口病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月達到第1個峰值，隨后第9—11月達到第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)性檢驗和白噪聲檢驗 對原始序列依次進行平穩(wěn)性檢驗和白噪聲檢驗，經(jīng)ADF檢驗發(fā)現(xiàn)，序列存在單位根(χ2=-0.018，P>0.05，見表2)，對原始序列的自相關圖、偏自相關圖進行分析，自相關系數(shù)超過2倍標準差(見圖2)，綜合考慮上述結果，認為原始序列為不平穩(wěn)序列。對序列進行季節(jié)性差分得到平穩(wěn)序列(χ2=-3.091，P<0.05)，隨后對該序列進行Ljung-Box檢驗，序列滿足非白噪聲要求(χ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ù)既往文獻經(jīng)驗[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)，對模型的殘差序列進行Ljung-Box檢驗，結果為白噪聲(P=0.510，見表3)，說明該模型數(shù)據(jù)提取完全，效果滿意。

表2 序列的平穩(wěn)性及白噪聲檢驗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模型參數(shù)估計和擬合優(yōu)度統(tǒng)計量結果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)絡

因為手足口病的數(shù)據(jù)周期為12個月，因此LSTM神經(jīng)網(wǎng)絡的窗口長度設置為12，即輸入層節(jié)點數(shù)為12，預測下個月手足口病發(fā)病率，輸出層節(jié)點數(shù)為1。

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

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

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

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

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

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

表5 隱藏層層數(shù)對模型預測性能的影響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神經(jīng)網(wǎng)絡的預測性能高于SARIMA模型，表明LSTM神經(jīng)網(wǎng)絡能更好地預測深圳市寶安區(qū)HFMD發(fā)病趨勢。

黑色虛線左側為兩種模型的擬合效果，右側為兩種模型的預測效果圖3 SARIMA和LSTM模型對2010年1月至2018年12月深圳市寶安區(qū)的HFMD月發(fā)病率擬合及預測效果的比較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年我國將手足口病納入丙類傳染病管理以來，手足口病的發(fā)病率一直位于我國丙類傳染病的首位[17]。手足口病是由不同腸道病毒感染所致，不同病原體之間存在差異，易在托幼機構及學校呈暴發(fā)流行，這類聚集性疫情的發(fā)生給國家和地區(qū)的衛(wèi)生保健系統(tǒng)造成較大的經(jīng)濟負擔[18]。本研究中，深圳市寶安區(qū)2014—2019年手足口病發(fā)病率顯著高于2009—2013年，可能原因跟該地區(qū)傳染病疫情監(jiān)測系統(tǒng)逐漸完善、外來人口數(shù)增多、人口流動性強等密切相關。因此，了解手足口病流行趨勢，構建恰當?shù)氖肿憧诓☆A測模型，可為相關部門制定手足口病防控策略提供科學依據(jù)。

本研究構建SARIMA模型時，先對原始序列進行平穩(wěn)性及白噪聲檢驗，判斷原始序列是否滿足平穩(wěn)非白噪聲，其次觀察序列的ACF、PACF圖確定模型的參數(shù)范圍，通過Python網(wǎng)格搜索自動擬合較佳參數(shù)，最后運用AIC最小原則，選取了SARIMA(0，0，2)(0，1，2)12為相對最優(yōu)模型。用該模型分別對深圳市寶安區(qū)的HFMD月發(fā)病率進行擬合及預測時，誤差指數(shù)較小，擬合及預測原始序列的趨勢性和周期性的效果較好，表明SARIMA模型預測深圳市寶安區(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模型的研究結果一致。

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

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