基于Python语言的ARIMA模型在天津市结核病发病率预测中的应用

张晓卉，姚婷婷，陈 阳，张甜甜，马 骏

(天津医科大学公共卫生学院流行病与卫生统计学系，天津 300041)

结核病是影响我国人民健康的重大传染病之一，发病率和病死率均位于我国传染病第二位[1-2]。世界卫生组织(WHO)2018年报告显示，结核病合并免疫系统受损(如感染人类免疫缺陷病毒)、营养不良、糖尿病以及吸烟等将大大提高其病死率[3]。因此，对结核病发病率进行预测，根据其未来预测趋势为相关部门制定防控措施，具有实际意义。目前，国内已有学者利用时间序列模型对结核病的发病趋势进行预测并取得了较好的效果[4-6]，但对天津市结核病月发病率预测时间序列研究鲜有报道。鉴于关于天津市结核病发病率水平的研究较少，本研究基于2004年1月—2016年12月天津市结核病月发病率数据，通过时间序列分析建立乘积季节差分自回归移动平均(autoregressive integrated moving average, ARIMA)模型，拟对天津市结核病月发病率进行预测，以期增强结核病防控预警。当前常见的统计分析软件为SAS、SPSS、R语言等，三者主要针对统计分析，对于大数据的挖掘开发可视化，用户普及性等不及Python语言，故本研究应用Python语言进行预测模型的拟合。

1 资料与方法

1.1 资料来源 天津市2004—2016结核病月发病率数据来源于传染病网络直报，数据获取平台为国家人口与健康科学数据共享平台公共卫生科学数据中心(http://www.ncmi.cn/info/69/1544)，资料可靠。

1.2 Python语言 Python语言作为当前最热门的编程语言之一，仅次于java语言和C语言[7]。Python语言不仅可用于统计分析，还被广泛的应用于数据爬取、人工智能等领域[8]，其语言简单、优美，且免费开源。强大的第三方库支持多种科学计算和统计分析[9]。在时间序列分析中，Python语言建模过程简单，图形直观。当处理的时间序列数据量较大时，Python语言可利用其第三方库pandas，从而规避循环，极大地节省程序运行时间，具有R语言等不具有的自身优势。

1.3 季节性差分自回归移动平均(SARIMA)模型 ARIMA模型将数据视为随时间变化的随机变量，根据序列的历史值预测将来值[10]。考虑到季节性，使用乘积SARIMA模型。模型常表述为SARIMA(p, d, q) (P, D, Q),其中AR是自回归，p为自回归项，MA为移动平均，q为移动平均项数，d为使时间序列平稳的差分次数，“P,D,Q”为相应带有季节性的参数。

1.3.1 建模步骤[11-14](1)平稳性检验：依据时间序列图、自相关图 (autocorrelation function, ACF)、偏自相关图 (partial autocorrelation function, PACF)及迪基-福勒检验(Dickey-Fuller Test)判断数据是否平稳，若不平稳需要进行平稳性处理。(2)平稳性处理：可采取对数并差分、季节差分、移动平均等方法使序列平稳。(3)白噪声检验：当BOX-Ljung统计量P≤0.05时，判断序列为非白噪声序列，可进行后续分析。(4)定阶：根据差分次数确定d、D，根据ACF和PACF确定p、q、P、Q，参考赤池信息准则(Akaike information criterion，AIC)及贝叶斯信息准则(Bayesian information criterion，BIC)确定最优模型。(5)参数估计及残差分析：确定模型参数并检验其显著性，对残差进行白噪声检验。若检验通过，则模型建模良好。

1.3.2 模型预测及评价 以天津市2004年1月—2015年12月结核病月发病率数据作为训练集拟合模型，以 2016年1—12月数据作为验证集验证模型拟合精度，预测2017年1月—2019年12月天津市结核病月发病率[15]。以误差均方 (mean square error, MSE)、平均绝对误差 (mean absolute error, MAE)衡量模型的预测精度和建模效果，其中MSE、MAE越小, 表示预测精度越好, 建模效果越好[16]。

1.4 统计学方法 建模平台采用基于64位Anaconda(Python3.7)，调用的模块主要有numpy、pandas、matplotlib以及statsmodels。

Python建模过程，(1)导入一系列包：import pandas as pd,import numpy as np，import statsmodels.api as sm，import matplotlib.pyplot as plt；(2)导入2004年1月—2015年12月数据并将其转换为时间类型数据：data=pd.read_csv('jiehebin.csv'),dateparse=lambda dates: pd.datetime.strptime(dates,'%Y-%m'),data=pd.read_csv('jiehebin.csv',parse_dates=['month'],index_col='month',date_parser=dateparse)；(3)差分使序列平稳：data_first_difference= data.diff(1)，data_seasonal_first_difference=data_first_difference.diff(12)；(4)模型拟合及预测：Mod=sm.tsa.statespace.SARIMAX(data,trend="n",order=(1,1,1),seasonal_order =(3,1,1,12))，results=mod.fit()，predict_ts=results.predict()。

2 结果

2.1 流行病学趋势 以2004年1月—2015年12月数据建模，将时间序列图分解为趋势性成分、季节性成分和随机成分三部分，结果显示2004年1月—2015年12月天津市结核病月发病率总体呈下降趋势，2005—2008年出现一个发病高峰，2009年后大幅度下降，随后趋于平稳。结核病发病率于每年3—8月份高发，冬季骤然降低，提示其具有季节性。见图1、2。

图2 天津市2004年1月—2015年12月结核病月发病率时间序列分解图

2.2 SARIMA模型建模结果

2.2.1 数据预处理 原始时间序列图经Dickey-Fuller Test，结果显示P值为0.81，原始序列不平稳，需要进行平稳性处理。对原始序列进行一次差分和一次季节差分后序列图平稳，见图3。ACF和PACF亦显示数据已平稳，见图4。Dickey-Fuller Test结果显示P值为0.000029，亦说明序列已平稳。经差分后的平稳序列BOX-Ljung统计量P<0.05，判断该序列为非白噪声序列。

图3 差分后序列时间序列图及滚动均值、滚动标准差图

图4 差分后序列的ACF、PACF

2.2.2 模型识别 由序列一阶差分和一阶季节差分初步确定d=1，D=1,初步确定模型形式为SARIMA(p,1,q)×(P,1,Q)12。观察差分后的ACF和PACF(见图4)，可见ACF呈截尾，PACF呈拖尾，提示季节性模型为SARIMA(3，1，1)12模型。观察SARIMA(1,1,1)×(3,1,1)12残差序列的ACF和PACF(见图5)，提示非季节模型为SARIMA(1，1，1)，故原始序列初步拟合为乘积混合效应模型SARIMA(1,1,1)×(3,1,1)12。为确保筛选出最优模型，采用从低阶到高阶逐个进行尝试的方法挑选模型参数[17-18]。初步纳入SARIMA(3,1,1)×(3,1,1)12、SARIMA(2,1,1)×(3,1,1)12、SARIMA(1,1,1)×(3,1,1)12、SARIMA(1,1,1)×(3,1,2)12、SARIMA(2,1,1)×(3,1,2)12、SARIMA(3,1,1)×(3,1,2)12六个模型进行试验，根据AIC、BIC准则选取其中最小值作为最优模型(见表1)，因此原始序列拟合为SARIMA(1,1,1)×(3,1,1)12，拟合后残差的ACF、PACF见图5。残差BOX-Ljung统计量P值为0.493,判断模型拟合后残差为白噪声序列。

表1 拟纳入模型AIC、BIC值

图5 SARIMA(1,1,1)×(3,1,1)12模型拟合后残差的ACF、PACF

2.2.3 参数估计及检验 模型的参数估计结果除ar.L1无统计学意义外，其他参数均有统计学意义。因此，除去ar.L1，将其他参数全部列入SARIMA(1,1,1)×(3,1,1)12模型。见表2。

表2 SARIMA(1,1,1)×(3,1,1)12模型的参数估计

2.2.4 模型拟合 将天津市2004年1月—2015年12月结核病月发病率数据作为训练集拟合SARIMA(1,1,1)×(3,1,1)12模型，其中MAE=0.306，MSE=0.224。见图6。

图6 2004年1月—2015年12月天津市结核病月发病率拟合结果

2.2.5 模型效果评价 将2016年1—12月结核病月发病率进行回代预测，实际发病例数均在预测发病例数的95%CI内，模型拟合良好，具有较好的预测性能，其中MAE=0.169，MSE=0.081，可对天津市结核病的发病数进行较准确的预测。见图7、表3。

图7 2016年1—12月天津市结核病月发病率预测结果

表3 2016年1—12月实际发病率与预测值比较(/10万)

2.2.6 模型预测 利用SARIMA(1,1,1)×(3,1,1)12对天津市2017年1月—2019年12月肺结核发病率进行预测，结果显示天津市结核病月发病率将总体呈现下降趋势，春季高发，冬季发病率降低，符合结核病发病规律，预测结果可供参考。见图8、表4。

图8 SARIMA对2017年1月—2019年12月天津市结核病月发病率预测结果

表4 SARIMA对2019年7—12月天津市结核病月发病的预测值

3 讨论

随着深度学习和数据挖掘技术的日渐成熟，Python语言风靡全球。在信息爆炸和多学科融合的大数据时代，Python语言作为一门通用语言在统计学应用中的地位也愈加重要，相较于R语言，Python语言在前期数据收集、处理、建模和运行速度等方面显示出卓越性能，因此本研究使用Python语言建立ARIMA时间序列预测模型。

目前，有多种模型可用于传染病发病率的预测，如GM(1,1)模型[19]、马尔可夫链预测模型[20]、ARIMA模型[13,21-22]等。其中ARIMA模型作为最经典的预测模型之一，可将时间序列分解为趋势性成分、季节性成分和随机干扰，并对噪声进行分析处理，预测精度较高，是传染病时间序列预测模型中最重要的手段[23]。胡晓媛等[24]应用SARIMA模型对全国肺结核月发病率进行预测，预测MAE值为0.416992。本研究MAE值为0.169,较之稍高。秘玉清等[25]应用SARIMA模型预测山东省结核病发病趋势，得出发病率将呈现周期性上升的结论。鉴于对天津市肺结核月发病率研究较少，本研究对天津市结核病月发病率的流行病学趋势作了描述性研究，并预测其未来发病趋势，可为天津市相关部门制定防控措施提供理论依据。

本研究基于2004年1月—2016年12月天津市结核病月发病率资料，将时间序列分解为趋势性、季节性和随机噪声三部分，分析结核病的发病趋势和季节性，最终确定SARIMA(1,1,1)×(3,1,1)12为天津市结核病月发病率的最终模型。首先，利用 2004年1月—2015年12月数据建立最优模型，结果显示SARIMA(1,1,1)×(3,1,1)12可较准确地拟合实际月发病率，残差为白噪声序列，说明建模良好。然后将2016年1—12月结核病月发病率进行回代预测，实际值均在预测值95%CI内，说明模型适用于对天津市未来结核病月发病率的预测。最后，将SARIMA(1,1,1)×(3,1,1)12模型应用于对2017年1月—2019年12月结核病月发病率的预测，可用预测值来估计未来结核病的流行强度，若2017年1月—2019年12月实际发病人数在预测发病人数的95%CI内, 表明当月的结核病疫情基本正常；若发现实际发病人数处于预测发病人数的95%CI外, 则提示结核病疫情有可能发生异常[15]。

本研究将2016年1—12月发病率进行回代预测时，出现8—10月份发病率相对误差和其他月份相差较大的情况，其原因：一方面，可能是由于ARIMA模型作为经典的线性模型在处理具有非线性特点的时间序列问题上表现出一定的局限性；另一方面，结核病的发生发展具有多因素性，未来考虑将除时间因素以外的其他混杂因素列入模型，以提高其预测精度。

时间序列具有混沌现象，存在内在随机性和不规则有序性[11]。因此，对时间序列的预测在尽可能抓取线性成分的基础上还应更多的关注非线性成分。而ARIMA模型作为一种线性模型在处理非线性成分的问题上具有一定的不足，容易使预测精度降低。目前，关于基于误差补偿思想和相空间重构思想的ARIMA混合模型已经崭露头角[26-27]。然而，对于ARIMA-SVM、ARIMA-LSTM等混合模型虽预测效果较单纯ARIMA模型较好，但解释性差。在ARIMA模型建模过程中，残差应通过白噪声检验才能判断建模是否合理，而对残差进行相空间重构则需要残差具有混沌性，而非随机序列，以上矛盾之处此类混合模型无法给出合理的解释。另外，关于将ARIMA模型和各类神经网络相结合的混合模型也层出不穷，研究结果显示其预测精度相较于单纯ARIMA模型高，但由于神经网络背后数学理论仍处于“黑箱子”状态，其混合模型预测的准确性和重复性仍有待探讨。神经网络模型的训练效果和预测精度随着数据量的增大和数据维度的增高而增大，基于时间序列数据短期预测的特点，将神经网络运用于时间序列中，其理论有待完善。ARIMA模型常用于短期预测，因此，其预测结果需要随着数据量的更新而不断更新，如何将ARIMA和其他各种预测模型如灰色模型、隐马尔可夫链等有效结合，通过充分提取时间序列的线性和非线性成分，控制随机误差，提高预测的准确度和精确度，是今后研究的方向。