ARIMA、ARIMAX、NGO-LSTM模型在山东省聊城市结核病发病数预测中的应用

孙明浩 段雨琪 郑良 于胜男 程传龙 左慧 陈鸣 李秀君

结核病是由结核分枝杆菌(Mycobacteriumtuberculosis,MTB)感染引起的一种传染性疾病,肺结核是最常见的类型[1]。根据《2022年全球结核病报告》,在30个结核病高负担国家中,中国估算结核病发病数位居第三位[2]。聊城市作为山东省结核病发病率较高的地区,防治形势较为严峻[3]。为有效防控结核病,世界卫生组织于2014年提出了终止结核病战略,目标是在2015—2035年期间将结核病死亡例数减少95%,并将发病例数减少90%[4]。

为了实现这一目标,准确预测结核病发病趋势显得十分重要。传染病预测模型主要有传统预测模型、基于机器学习的预测模型、动力学模型、空间模型和网络模型等[5-9]。传统预测模型有差分自回归移动平均(autoregressive integrated moving average,ARIMA)模型[10]、具有外生回归变量的差分自回归移动平均(autoregressive integrated moving average with exogenous regressors,ARIMAX)模型[11]、指数平滑模型[12]等;基于机器学习的预测模型有BP人工神经网络模型[13]、支持向量机(support vector machines,SVM)[14]、长期短期记忆(long short-term memory,LSTM)神经网络[15]等。目前,ARIMA和ARIMAX模型已被应用于结核病[11, 16]、手足口病[17-18]和新型冠状病毒感染[19]等疾病的预测中。同时,以往研究表明,温度、相对湿度、气压、风速、降水等气象因素均可能会影响结核病的发生[20-22]。一项中国东部和中国宁波市的时间序列研究也表明,ARIMA模型在纳入气象因素后预测性能有所改善[11, 23],可见纳入气象因素的ARIMAX模型可能会提高ARIMA模型的预测性能。另外,基于机器学习的预测模型在处理时间序列中的非线性关系方面具有良好的性能[15],并且可以结合优化算法得到较好的模型预测效果[24]。

基于此,本研究尝试比较ARIMA模型、考虑外部气象变量的ARIMAX模型和结合北方苍鹰优化(Northern goshawk optimization,NGO)算法的LSTM神经网络模型对聊城市结核病患者例数的预测效果,确定预测聊城市结核病流行趋势的最优预测模型,为有关部门针对结核病的预防及控制决策提供科学依据。

资料和方法

一、研究地区和资料收集

聊城市(如图1所示)位于山东省西部,介于东经115°16′—116°32′和北纬35°47′—37°02′之间,总面积8628 km2,截至2021年末,聊城市常住人口为592.79万人。

图1 聊城市在山东省的位置分布图[审图号GS(2022)1873号]

本研究患者数据来源于“中国疾病预防控制信息系统”子系统“结核病系统”。提取2011年1月至2018年12月的结核病月发病数作为研究对象。

气象因素包括月平均气温(MAT,℃)、月平均最高气温(MAHT,℃)、月平均最低气温(MALT,℃)、月平均相对湿度(MAH,%)、月平均气压(MAP,hPa)和月平均风速(MAS,m/s),数据来自国家气象科学数据中心(https://data.cma.cn/)。

二、研究方法

2. ARIMAX模型:选择最优ARIMA模型后,考虑纳入气象因子。计算上述最优ARIMA模型残差序列与环境因子最优ARIMA模型残差序列之间的互相关函数(cross-correlation function,CCF),识别显著的时间滞后,为纳入的气象因子选择滞后期。基于先前研究的结果,以及生物学和流行病学知识,最大滞后期选择为12个月[11]。

3. NGO-LSTM神经网络模型:对于基于机器学习的预测模型,使用NGO-LSTM模型对结核病发病数进行拟合预测[16, 24]。LSTM模型最早由Hochreiter和Schmidhuber[26]于1997年提出,LSTM单元由输入门、遗忘门和输出门组成。LSTM选择性地允许信息通过门结构,以更新或保留历史信息。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算法包括勘探与开发两个阶段,分别对空间进行全局搜索和局部搜索[24]。NGO算法两阶段更新所有种群成员后,确定了种群成员的新值、目标函数和最佳建议解决方案,随后算法进入下一次迭代,直到算法到达最后一次迭代,引入算法迭代过程中得到的最优解作为给定优化问题的准最优解[24]。本研究利用NGO算法优化LSTM算法的3个参数:隐藏层节点数、初始学习率和L2正则化系数。以2011年1月至2017年12月的结核病发病数据构建NGO-LSTM模型,预测2018年1—12月的结核病发病数。设置LSTM模型最大迭代次数为2000次,第850次迭代开始调整学习率,学习率调整因子为0.2,优化器使用ADAM,设置Look_Back=12。NGO种群数量和迭代次数分别设置为10和30。使用Matlab 2021a软件完成本研究3种模型的拟合和预测。

5.敏感性分析:ARIMA(p,d,q)(P,D,Q)s模型中参数p、d、q和D均根据数据特征确定,因此改变参数P和Q进行敏感性分析。同时也以次优ARIMA模型[ARIMA(3,1,3)×(1,1,1)12模型]为基模型通过相同方法进行ARIMAX分析。在NGO-LSTM模型中,本研究改变LSTM模型最大迭代次数、调整学习率时的迭代次数、学习率调整因子和优化器及NGO种群数量和迭代次数进行敏感性分析。

结 果

一、描述性统计

表1列出了聊城市2011—2018年结核病月发病数和月平均气象因子的描述性统计数据。如图2所示,月发病数呈季节性波动,高峰期主要发生在3、6、12月,低谷期多见于1、2月。

表1 2011—2018年山东省聊城市结核病月发病数和月平均气象因子

注 MAP:月平均气压(hPa);MAS:月平均风速(m/s);MAT:月平均温度(℃);MAH:月平均相对湿度(%);MALT:月平均最低温度(℃);MAHT:月平均最高温度(℃)图2 2011—2018年山东省聊城市结核病月发病数和气象因子时间序列

二、ARIMA模型

ARIMA模型对数据的基本要求是数据平稳,经过一阶差分和一阶季节差分(图3),时间序列平稳,因此参数d和D均为1。用自相关函数(autocorrelation function,ACF)和偏自相关函数(partial autocorrelation function,PACF)图来确定ARIMA模型参数p、q、P和Q的值。如图4和图5所示,自相关图和偏自相关图均是三阶截尾,p=3,q=3,参数P、Q的取值判定较为困难,因P、Q取值超过2阶的情况比较少见,故可在0、1、2中进行筛选,AIC值最小的模型作为最优模型。如表2所示,最优模型为ARIMA(3,1,3)×(0,1,1)12,AIC=751.211,Ljung-Box检验P值为0.881。

表2 候选ARIMA模型的比较

图3 一阶差分和一阶季节差分之后的时间序列曲线图

图4 结核病患者时间序列差分后自相关函数图

三、ARIMAX模型

本研究计算了结核病发病数和不同滞后时间的环境因素之间的CCF系数。表3显示,结核病发病数与MAT(滞后3个月,滞后5个月)、MAHT(滞后6个月)、MALT(滞后2个月)、MAH(滞后1个月)、MAP(滞后2个月)呈负相关,与MAS(滞后5个月)呈正相关。

表3 不同时滞的气象因素残差序列与结核病发病残差序列之间的互相关函数系数

如表4所示,具有MAH(滞后1个月)的ARIMA(3,1,3)×(0,1,1)12模型是最优的单变量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模型拟合和预测MAPE、MAE、RMSE均最小,在所有ARIMA和ARIMAX模型中拟合和预测效果最好,具有MAH(滞后1个月)的ARIMA(3,1,3)×(0,1,1)12模型次之。

表4 最优ARIMA模型、ARIMAX模型拟合预测效果比较

四、NGO-LSTM模型

以2011年1月至2017年12月的结核病病例数据构建NGO-LSTM模型,预测2018年结核病病例数据。NGO-LSTM模型具有较好的拟合和预测效果,MAPE(%)分别为5.510和5.820,MAE分别为13.045和13.119,RMSE分别为16.026和16.297。

五、模型比较

不同模型的预测结果如表5和图6所示。ARIMA模型、考虑月平均相对湿度(滞后1个月)与最低温度(滞后2个月)的多变量ARIMAX模型和NGO-LSTM模型对2018年结核病发病数预测的MAPE分别为9.293%、8.419%、5.820%,MAE分别为19.282、16.997、13.119,RMSE分别为23.773、22.191、16.297。NGO-LSTM神经网络模型的预测效果最好,MAPE最小,具有MAH(滞后1个月)和MALT(滞后2个月)的ARIMA(3,1,3)×(0,1,1)12模型预测效果次之。

表5 3种模型预测效果的比较

注 ARIMA:差分自回归移动平均模型;MAH:月平均相对湿度;MALT:月平均最低温度;lag1～lag2:滞后1～2个月;NGO-LSTM:结合北方苍鹰优化算法的长期短期记忆神经网络模型图6 3种模型预测值与实际值比较

六、敏感性分析

如表6所示,ARIMA模型中参数和的改变和ARIMAX模型中ARIMA部分的改变均对预测效果影响不大,结果较稳健。如表7所示,NGO-LSTM模型参数的改变对预测效果影响不大,结果较稳健。

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

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

讨 论

本研究通过构建ARIMA、ARIMAX和NGO-LSTM模型,比较3种模型对聊城市结核病发病数的预测效果。NGO-LSTM模型得到较好的预测效果,这可能是因为NGO优化算法根据数据特征优化了隐藏层节点数、初始学习率和L2正则化系数这3个超参数。隐藏层节点数主要影响模型的复杂度和容量,隐藏层数量越多,模型表达能力越强[27]。初始学习率决定模型优化器梯度下降时权重和偏差更新幅度,过大可能会导致训练不稳定或梯度爆炸的问题,过小可能会使得模型收敛速度变慢[28]。L2正则化可以防止过拟合,L2正则化系数决定了正则化项在损失函数中所占的比重,如果L2正则化系数过大,模型可能会导致欠拟合,如果过小,则无法有效约束模型,可能会导致模型过拟合[29]。并且NGO算法能够为优化问题提供理想的准最优解,其在优化方面的表现也明显优于粒子群优化算法和遗传算法等竞争算法[24]。综上,NGO-LSTM能够得到较理想的预测效果。

ARIMA和ARIMAX模型对结核病发病趋势的预测效果与以往研究相似[11]。NGO-LSTM模型对发病趋势的预测效果优于ARIMA和ARIMAX模型,这可能是因为ARIMA和ARIMAX模型能分析传染病序列的线性部分,但传染病数据的非线性部分可能不是白噪声, ARIMA和ARIMAX模型可能无法充分捕捉信息,而LSTM模型是一种深度学习应用程序,旨在学习时间模式、捕获非线性依赖关系并存储有用记忆,因此,NGO-LSTM模型可能会产生更好的结果[17, 26]。本研究利用ARIMA、ARIMAX和NGO-LSTM模型对聊城市结核病发病数进行预测,为其他地区类似研究提供了分析框架和思路,有助于地方政府和相关防控部门科学评估和预测结核病的流行动态,提前制定相应的防控措施,减少疾病的流行和传播。

本研究也存在一些局限性:数据来源于“中国疾病预防控制信息系统”子系统“结核病系统”,该系统为结核病登记数据,即在结核病定点医疗机构治疗管理的患者,可能与实际发病数据存在差异,这种差异体现在两个方面,一方面是报告的滞后性,另一方面是季节性人群行为(如春节、入学入职体检)等造成的影响;本研究只应用了ARIMA、ARIMAX和NGO-LSTM模型,未来可纳入其他模型并应用其他优化算法进行比较。

综上所述,NGO-LSTM模型对聊城市结核病月发病数的预测效果最好,可被用来预测结核病发病趋势,为有关部门针对结核病的预防及控制决策提供科学依据。

利益冲突所有作者均声明不存在利益冲突

作者贡献孙明浩:设计实施研究、撰写及修改文章;段雨琪、郑良、于胜男、程传龙和左慧:修改文章;陈鸣和李秀君:设计研究、对文章的知识性内容作批评性审阅、指导