ARIMA、ARIMAX、NGO-LSTM模型在山東省聊城市結(jié)核病發(fā)病數(shù)預(yù)測中的應(yīng)用

孫明浩 段雨琪 鄭良 于勝男 程傳龍 左慧 陳鳴 李秀君

結(jié)核病是由結(jié)核分枝桿菌(Mycobacteriumtuberculosis,MTB)感染引起的一種傳染性疾病,肺結(jié)核是最常見的類型[1]。根據(jù)《2022年全球結(jié)核病報告》,在30個結(jié)核病高負(fù)擔(dān)國家中,中國估算結(jié)核病發(fā)病數(shù)位居第三位[2]。聊城市作為山東省結(jié)核病發(fā)病率較高的地區(qū),防治形勢較為嚴(yán)峻[3]。為有效防控結(jié)核病,世界衛(wèi)生組織于2014年提出了終止結(jié)核病戰(zhàn)略,目標(biāo)是在2015—2035年期間將結(jié)核病死亡例數(shù)減少95%,并將發(fā)病例數(shù)減少90%[4]。

為了實現(xiàn)這一目標(biāo),準(zhǔn)確預(yù)測結(jié)核病發(fā)病趨勢顯得十分重要。傳染病預(yù)測模型主要有傳統(tǒng)預(yù)測模型、基于機(jī)器學(xué)習(xí)的預(yù)測模型、動力學(xué)模型、空間模型和網(wǎng)絡(luò)模型等[5-9]。傳統(tǒng)預(yù)測模型有差分自回歸移動平均(autoregressive integrated moving average,ARIMA)模型[10]、具有外生回歸變量的差分自回歸移動平均(autoregressive integrated moving average with exogenous regressors,ARIMAX)模型[11]、指數(shù)平滑模型[12]等;基于機(jī)器學(xué)習(xí)的預(yù)測模型有BP人工神經(jīng)網(wǎng)絡(luò)模型[13]、支持向量機(jī)(support vector machines,SVM)[14]、長期短期記憶(long short-term memory,LSTM)神經(jīng)網(wǎng)絡(luò)[15]等。目前,ARIMA和ARIMAX模型已被應(yīng)用于結(jié)核病[11, 16]、手足口病[17-18]和新型冠狀病毒感染[19]等疾病的預(yù)測中。同時,以往研究表明,溫度、相對濕度、氣壓、風(fēng)速、降水等氣象因素均可能會影響結(jié)核病的發(fā)生[20-22]。一項中國東部和中國寧波市的時間序列研究也表明,ARIMA模型在納入氣象因素后預(yù)測性能有所改善[11, 23],可見納入氣象因素的ARIMAX模型可能會提高ARIMA模型的預(yù)測性能。另外,基于機(jī)器學(xué)習(xí)的預(yù)測模型在處理時間序列中的非線性關(guān)系方面具有良好的性能[15],并且可以結(jié)合優(yōu)化算法得到較好的模型預(yù)測效果[24]。

基于此,本研究嘗試比較ARIMA模型、考慮外部氣象變量的ARIMAX模型和結(jié)合北方蒼鷹優(yōu)化(Northern goshawk optimization,NGO)算法的LSTM神經(jīng)網(wǎng)絡(luò)模型對聊城市結(jié)核病患者例數(shù)的預(yù)測效果,確定預(yù)測聊城市結(jié)核病流行趨勢的最優(yōu)預(yù)測模型,為有關(guān)部門針對結(jié)核病的預(yù)防及控制決策提供科學(xué)依據(jù)。

資料和方法

一、研究地區(qū)和資料收集

聊城市(如圖1所示)位于山東省西部,介于東經(jīng)115°16′—116°32′和北緯35°47′—37°02′之間,總面積8628 km2,截至2021年末,聊城市常住人口為592.79萬人。

圖1 聊城市在山東省的位置分布圖[審圖號GS(2022)1873號]

本研究患者數(shù)據(jù)來源于“中國疾病預(yù)防控制信息系統(tǒng)”子系統(tǒng)“結(jié)核病系統(tǒng)”。提取2011年1月至2018年12月的結(jié)核病月發(fā)病數(shù)作為研究對象。

氣象因素包括月平均氣溫(MAT,℃)、月平均最高氣溫(MAHT,℃)、月平均最低氣溫(MALT,℃)、月平均相對濕度(MAH,%)、月平均氣壓(MAP,hPa)和月平均風(fēng)速(MAS,m/s),數(shù)據(jù)來自國家氣象科學(xué)數(shù)據(jù)中心(https://data.cma.cn/)。

二、研究方法

2. ARIMAX模型:選擇最優(yōu)ARIMA模型后,考慮納入氣象因子。計算上述最優(yōu)ARIMA模型殘差序列與環(huán)境因子最優(yōu)ARIMA模型殘差序列之間的互相關(guān)函數(shù)(cross-correlation function,CCF),識別顯著的時間滯后,為納入的氣象因子選擇滯后期?；谙惹把芯康慕Y(jié)果,以及生物學(xué)和流行病學(xué)知識,最大滯后期選擇為12個月[11]。

3. NGO-LSTM神經(jīng)網(wǎng)絡(luò)模型:對于基于機(jī)器學(xué)習(xí)的預(yù)測模型,使用NGO-LSTM模型對結(jié)核病發(fā)病數(shù)進(jìn)行擬合預(yù)測[16, 24]。LSTM模型最早由Hochreiter和Schmidhuber[26]于1997年提出,LSTM單元由輸入門、遺忘門和輸出門組成。LSTM選擇性地允許信息通過門結(jié)構(gòu),以更新或保留歷史信息。LSTM可以表示為:

ft=σ(Wf·[ht-1,xt]+bf)#(1)

it=σ(Wi·[ht-1,xt]+bi)#(2)

ht=ot·tanh(Ct)#(5)

NGO算法于2021年由Mohammad等[24]提出。NGO算法包括勘探與開發(fā)兩個階段,分別對空間進(jìn)行全局搜索和局部搜索[24]。NGO算法兩階段更新所有種群成員后,確定了種群成員的新值、目標(biāo)函數(shù)和最佳建議解決方案,隨后算法進(jìn)入下一次迭代,直到算法到達(dá)最后一次迭代,引入算法迭代過程中得到的最優(yōu)解作為給定優(yōu)化問題的準(zhǔn)最優(yōu)解[24]。本研究利用NGO算法優(yōu)化LSTM算法的3個參數(shù):隱藏層節(jié)點數(shù)、初始學(xué)習(xí)率和L2正則化系數(shù)。以2011年1月至2017年12月的結(jié)核病發(fā)病數(shù)據(jù)構(gòu)建NGO-LSTM模型,預(yù)測2018年1—12月的結(jié)核病發(fā)病數(shù)。設(shè)置LSTM模型最大迭代次數(shù)為2000次,第850次迭代開始調(diào)整學(xué)習(xí)率,學(xué)習(xí)率調(diào)整因子為0.2,優(yōu)化器使用ADAM,設(shè)置Look_Back=12。NGO種群數(shù)量和迭代次數(shù)分別設(shè)置為10和30。使用Matlab 2021a軟件完成本研究3種模型的擬合和預(yù)測。

5.敏感性分析:ARIMA(p,d,q)(P,D,Q)s模型中參數(shù)p、d、q和D均根據(jù)數(shù)據(jù)特征確定,因此改變參數(shù)P和Q進(jìn)行敏感性分析。同時也以次優(yōu)ARIMA模型[ARIMA(3,1,3)×(1,1,1)12模型]為基模型通過相同方法進(jìn)行ARIMAX分析。在NGO-LSTM模型中,本研究改變LSTM模型最大迭代次數(shù)、調(diào)整學(xué)習(xí)率時的迭代次數(shù)、學(xué)習(xí)率調(diào)整因子和優(yōu)化器及NGO種群數(shù)量和迭代次數(shù)進(jìn)行敏感性分析。

結(jié) 果

一、描述性統(tǒng)計

表1列出了聊城市2011—2018年結(jié)核病月發(fā)病數(shù)和月平均氣象因子的描述性統(tǒng)計數(shù)據(jù)。如圖2所示,月發(fā)病數(shù)呈季節(jié)性波動,高峰期主要發(fā)生在3、6、12月,低谷期多見于1、2月。

表1 2011—2018年山東省聊城市結(jié)核病月發(fā)病數(shù)和月平均氣象因子

注 MAP:月平均氣壓(hPa);MAS:月平均風(fēng)速(m/s);MAT:月平均溫度(℃);MAH:月平均相對濕度(%);MALT:月平均最低溫度(℃);MAHT:月平均最高溫度(℃)圖2 2011—2018年山東省聊城市結(jié)核病月發(fā)病數(shù)和氣象因子時間序列

二、ARIMA模型

ARIMA模型對數(shù)據(jù)的基本要求是數(shù)據(jù)平穩(wěn),經(jīng)過一階差分和一階季節(jié)差分(圖3),時間序列平穩(wěn),因此參數(shù)d和D均為1。用自相關(guān)函數(shù)(autocorrelation function,ACF)和偏自相關(guān)函數(shù)(partial autocorrelation function,PACF)圖來確定ARIMA模型參數(shù)p、q、P和Q的值。如圖4和圖5所示,自相關(guān)圖和偏自相關(guān)圖均是三階截尾,p=3,q=3,參數(shù)P、Q的取值判定較為困難,因P、Q取值超過2階的情況比較少見,故可在0、1、2中進(jìn)行篩選,AIC值最小的模型作為最優(yōu)模型。如表2所示,最優(yōu)模型為ARIMA(3,1,3)×(0,1,1)12,AIC=751.211,Ljung-Box檢驗P值為0.881。

表2 候選ARIMA模型的比較

圖3 一階差分和一階季節(jié)差分之后的時間序列曲線圖

圖4 結(jié)核病患者時間序列差分后自相關(guān)函數(shù)圖

三、ARIMAX模型

本研究計算了結(jié)核病發(fā)病數(shù)和不同滯后時間的環(huán)境因素之間的CCF系數(shù)。表3顯示,結(jié)核病發(fā)病數(shù)與MAT(滯后3個月,滯后5個月)、MAHT(滯后6個月)、MALT(滯后2個月)、MAH(滯后1個月)、MAP(滯后2個月)呈負(fù)相關(guān),與MAS(滯后5個月)呈正相關(guān)。

表3 不同時滯的氣象因素殘差序列與結(jié)核病發(fā)病殘差序列之間的互相關(guān)函數(shù)系數(shù)

如表4所示,具有MAH(滯后1個月)的ARIMA(3,1,3)×(0,1,1)12模型是最優(yōu)的單變量ARIMAX模型,具有MALT(滯后2個月)的ARIMA(3,1,3)×(0,1,1)12模型次之,其MAPE、MAE、RMSE均小于ARIMA模型。將MALT(滯后2個月)和MAH(滯后1個月)作為外生回歸變量納入多變量ARIMAX模型。具有MALT(滯后2個月)和MAH(滯后1個月)的ARIMA(3,1,3)×(0,1,1)12模型擬合和預(yù)測MAPE、MAE、RMSE均最小,在所有ARIMA和ARIMAX模型中擬合和預(yù)測效果最好,具有MAH(滯后1個月)的ARIMA(3,1,3)×(0,1,1)12模型次之。

表4 最優(yōu)ARIMA模型、ARIMAX模型擬合預(yù)測效果比較

四、NGO-LSTM模型

以2011年1月至2017年12月的結(jié)核病病例數(shù)據(jù)構(gòu)建NGO-LSTM模型,預(yù)測2018年結(jié)核病病例數(shù)據(jù)。NGO-LSTM模型具有較好的擬合和預(yù)測效果,MAPE(%)分別為5.510和5.820,MAE分別為13.045和13.119,RMSE分別為16.026和16.297。

五、模型比較

不同模型的預(yù)測結(jié)果如表5和圖6所示。ARIMA模型、考慮月平均相對濕度(滯后1個月)與最低溫度(滯后2個月)的多變量ARIMAX模型和NGO-LSTM模型對2018年結(jié)核病發(fā)病數(shù)預(yù)測的MAPE分別為9.293%、8.419%、5.820%,MAE分別為19.282、16.997、13.119,RMSE分別為23.773、22.191、16.297。NGO-LSTM神經(jīng)網(wǎng)絡(luò)模型的預(yù)測效果最好,MAPE最小,具有MAH(滯后1個月)和MALT(滯后2個月)的ARIMA(3,1,3)×(0,1,1)12模型預(yù)測效果次之。

表5 3種模型預(yù)測效果的比較

注 ARIMA:差分自回歸移動平均模型;MAH:月平均相對濕度;MALT:月平均最低溫度;lag1～lag2:滯后1～2個月;NGO-LSTM:結(jié)合北方蒼鷹優(yōu)化算法的長期短期記憶神經(jīng)網(wǎng)絡(luò)模型圖6 3種模型預(yù)測值與實際值比較

六、敏感性分析

如表6所示,ARIMA模型中參數(shù)和的改變和ARIMAX模型中ARIMA部分的改變均對預(yù)測效果影響不大,結(jié)果較穩(wěn)健。如表7所示,NGO-LSTM模型參數(shù)的改變對預(yù)測效果影響不大,結(jié)果較穩(wěn)健。

表6 ARIMA和ARIMAX模型的敏感性分析

表7 NGO-LSTM模型的敏感性分析

討 論

本研究通過構(gòu)建ARIMA、ARIMAX和NGO-LSTM模型,比較3種模型對聊城市結(jié)核病發(fā)病數(shù)的預(yù)測效果。NGO-LSTM模型得到較好的預(yù)測效果,這可能是因為NGO優(yōu)化算法根據(jù)數(shù)據(jù)特征優(yōu)化了隱藏層節(jié)點數(shù)、初始學(xué)習(xí)率和L2正則化系數(shù)這3個超參數(shù)。隱藏層節(jié)點數(shù)主要影響模型的復(fù)雜度和容量,隱藏層數(shù)量越多,模型表達(dá)能力越強(qiáng)[27]。初始學(xué)習(xí)率決定模型優(yōu)化器梯度下降時權(quán)重和偏差更新幅度,過大可能會導(dǎo)致訓(xùn)練不穩(wěn)定或梯度爆炸的問題,過小可能會使得模型收斂速度變慢[28]。L2正則化可以防止過擬合,L2正則化系數(shù)決定了正則化項在損失函數(shù)中所占的比重,如果L2正則化系數(shù)過大,模型可能會導(dǎo)致欠擬合,如果過小,則無法有效約束模型,可能會導(dǎo)致模型過擬合[29]。并且NGO算法能夠為優(yōu)化問題提供理想的準(zhǔn)最優(yōu)解,其在優(yōu)化方面的表現(xiàn)也明顯優(yōu)于粒子群優(yōu)化算法和遺傳算法等競爭算法[24]。綜上,NGO-LSTM能夠得到較理想的預(yù)測效果。

ARIMA和ARIMAX模型對結(jié)核病發(fā)病趨勢的預(yù)測效果與以往研究相似[11]。NGO-LSTM模型對發(fā)病趨勢的預(yù)測效果優(yōu)于ARIMA和ARIMAX模型,這可能是因為ARIMA和ARIMAX模型能分析傳染病序列的線性部分,但傳染病數(shù)據(jù)的非線性部分可能不是白噪聲, ARIMA和ARIMAX模型可能無法充分捕捉信息,而LSTM模型是一種深度學(xué)習(xí)應(yīng)用程序,旨在學(xué)習(xí)時間模式、捕獲非線性依賴關(guān)系并存儲有用記憶,因此,NGO-LSTM模型可能會產(chǎn)生更好的結(jié)果[17, 26]。本研究利用ARIMA、ARIMAX和NGO-LSTM模型對聊城市結(jié)核病發(fā)病數(shù)進(jìn)行預(yù)測,為其他地區(qū)類似研究提供了分析框架和思路,有助于地方政府和相關(guān)防控部門科學(xué)評估和預(yù)測結(jié)核病的流行動態(tài),提前制定相應(yīng)的防控措施,減少疾病的流行和傳播。

本研究也存在一些局限性:數(shù)據(jù)來源于“中國疾病預(yù)防控制信息系統(tǒng)”子系統(tǒng)“結(jié)核病系統(tǒng)”,該系統(tǒng)為結(jié)核病登記數(shù)據(jù),即在結(jié)核病定點醫(yī)療機(jī)構(gòu)治療管理的患者,可能與實際發(fā)病數(shù)據(jù)存在差異,這種差異體現(xiàn)在兩個方面,一方面是報告的滯后性,另一方面是季節(jié)性人群行為(如春節(jié)、入學(xué)入職體檢)等造成的影響;本研究只應(yīng)用了ARIMA、ARIMAX和NGO-LSTM模型,未來可納入其他模型并應(yīng)用其他優(yōu)化算法進(jìn)行比較。

綜上所述,NGO-LSTM模型對聊城市結(jié)核病月發(fā)病數(shù)的預(yù)測效果最好,可被用來預(yù)測結(jié)核病發(fā)病趨勢,為有關(guān)部門針對結(jié)核病的預(yù)防及控制決策提供科學(xué)依據(jù)。

利益沖突所有作者均聲明不存在利益沖突

作者貢獻(xiàn)孫明浩:設(shè)計實施研究、撰寫及修改文章;段雨琪、鄭良、于勝男、程傳龍和左慧:修改文章;陳鳴和李秀君:設(shè)計研究、對文章的知識性內(nèi)容作批評性審閱、指導(dǎo)