ARIMA乘积季节模型在青州市布鲁氏菌病发病预测中的应用

刘 杰 武钦发 李伟国 肖宇飞 毛 倩 石福艳 王素珍

1 山东省青州市地方病防治研究所，262500 山东 青州； 2 山东省青州市疾病预防控制中心，262500 山东 青州；3 潍坊医学院卫生统计学教研室，261053 山东 潍坊

布鲁氏菌病是由布鲁氏菌属细菌入侵机体引起人和牛、羊、猪等动物共患的传染病[1]。人布鲁氏菌病被《中华人民共和国传染病防治法》列入乙类法定报告传染病管理。青州市首次在1965年发现布鲁氏菌病病人，其后布鲁氏菌病在青州市的流行基本处于稳定状态，但是2011年以来,由于养殖业的不断发展、牲畜交易市场开放、相应的防治措施落实不到位，布鲁氏菌病疫情有不断上升趋势。为了更好地掌握青州市布鲁氏菌病流行趋势，本文根据青州市2011—2017年布鲁氏菌病月发病数建立自回归求和移动平均(ARIMA)乘积季节模型，并应用该模型预测青州市2018年各月份布鲁氏菌病发病数，用实际月发病数进行回代拟合，检验模型的预测效果,以探讨ARIMA 模型预测布鲁氏菌病月发病数的可行性。

1 资料与方法

1.1 资料来源

青州市2011—2017年布鲁氏菌病疫情资料来源于中国疾病预防控制系统，人口学资料来源于青州市统计年鉴和2011—2017年青州市国民经济和社会发展公报。

1.2 ARIMA模型简介

1.2.1定义

ARIMA模型全称为自回归移动平均模型，记作ARIMA(p,d,q)，由西方学者George Edward Pelham Box和 Gwilym Meirion Jenkins于20世纪70年代提出，是一种著名的时间序列统计模型[2]。其含义是指将非平稳时间序列转化为平稳时间序列，然后将因变量仅对它的滞后值以及随机误差项的现值和滞后值进行回归所建立的模型。

1.2.2建模步骤

第1步，时间序列平稳化。首先绘制时间序列的散点图或折线图，粗略判断其平稳性。对于非平稳的时间序列，如果存在异方差，则应对序列进行对数转换；如果时间序列存在一定变化趋势，则应对序列进行差分处理，差分包括季节性差分和非季节性差分，一般取d与D值为1或2，使其转为平稳时间序列[2-3]。转换后的时间序列是否平稳可通过单位根检验(augmented dickey-fuller test, ADF)来判定。

第2步，模型识别与定阶。模型的识别就是根据时间序列的识别规则判断模型适用于哪种过程，建立拟合效果最优的模型。用到的工具为自相关函数ACF和偏自相关函数PACF以及它们各自的相关图。通过观察分析自相关函数和偏自相关函数的拖尾性和结尾性判断时间序列适用哪种模型[3-4]。对于季节性ARIMA模型，先识别非季节成分，确定p、q值，再识别季节成分，P、Q阶数较难确定，一般来说超过2 阶的情况很少见。可将阶数从低阶到高阶进行不同的组合，选择出几个粗模型[3,5]。

第3步，参数估计与诊断。在得到几个备选模型后，需进行参数估计，检验是否有统计学意义。参数估计的方法有极大似然估计、最小二乘法估计、矩估计等统计方法[2]。经过参数估计得到几个备选模型后，用Ljung-Box 统计量对模型进行适应性检验，检查该模型时间序列所蕴含的信息提取是否充分，若P>0.05则该模型残差是白噪声。同时根据AIC 和BIC 准则判定模型的拟合优度。经过反复试验比较，将模型系数有统计学意义、残差序列是白噪声、AIC 和BIC值较小的模型定为预测模型[5-6]。

第4步，预测及应用。利用选择出的最优模型进行预测，将预测值与实际值进行比较，评价模型的预测效果。通过计算平均绝对百分比误差(MAPE)、均方根误差(RMSE)等指标来评估模型预测的准确性和精度[7-8]。

1.3 统计学方法

本文采用Excel 2013 对2011—2017年青州市布鲁氏菌病月发病数据进行整理，采用SPSS 25.0软件构建青州市布鲁氏菌病的ARIMA 模型。α=0.05。

2 结果

2.1 序列平稳化

本研究以青州市2011—2017年布鲁氏菌病月发病数为基础数据，用2018年1—12月份青州市布鲁氏菌病实际月发病数对构建的ARIMA模型进行验证。首先绘制时间序列图(图1)，图1显示青州市布鲁氏菌病发病数总体呈现上升趋势，且具备季节性和周期性特点，发病集中在4—7月份，5月份为发病高峰，提示该时间序列是不平稳的时间序列，需要对该时间序列平稳化处理。由于存在周期性，对该序列同时进行一般性差分和季节性差分，周期为12，即d和D的取值为1。经差分后的时间序列图见图2，此时序列已基本平稳。

2.2 模型识别和参数估计

对获得的平稳时间序列绘制自相关函数图ACF和偏自相关函数图PACF(见图3)。由ACF图和PACF图可见，从一阶开始后均视为拖尾，因此建立ARIMA(p,1,q)(P,1,Q)12模型，其中p和q取值为0或1；根据经验，P和Q超过2阶的情况很少，采取从低阶到高阶依次组合的方法逐一进行尝试，应用Box-Ljung 方法对残差序列进行白噪声检验剔除非白噪声模型 (P<0.05) ，初步筛选出4 种平稳的白噪声模型ARIMA 模型。这4种模型的标准化BIC值和Box-Ljung 检验结果见表1。经过对比，这4种模型中ARIMA(0,1,1)(1,1,0)12模型BIC值相对较小，固定R平方相对较大并且该模型参数均具有统计学意义(P<0.01),该模型参数估计见表2。确定模型ARIMA(0,1,1)(1,1,0)12为最优模型。绘制ARIMA(0,1,1)(1,1,0)12模型的残差序列 ACF 和PACF 图(图4)对模型进行诊断，本模型残差序列 ACF 和 PACF 图显示残差基本都落在 95%区间内(P>0.05)。表明该模型的信息提取充分，残差为白噪声。初步认为该模型包含原始时间序列的所有特征，并且时间序列之间不存在相关性，可以用来预测青州市2018年各个月份人布鲁氏菌病月发病报告数。

图1 青州市2011—2017年布鲁氏菌病月发病人数序列图

图2 经过差分后青州市2011—2017年布鲁氏菌病月发病人数序列

图3 经过差分转换后的青州市布鲁氏菌病月发病数自相关图和偏自相关图

表1 4种备选ARIMA模型的BIC值和Ljung-Box Q结果

表2 ARIMA(0,1,1)(1,1,0)12模型的参数估计

图4 残差序列的自相关(ACF)和偏自相关(PACF)图

2.3 模型预测及应用

采用青州市2011年1月至2017年12月布鲁氏菌病月发病人数回代入ARIMA(0,1,1)(1,1,0)12模型进行拟合，绘制拟合曲线，见图5。结果显示拟合值与观测值总体趋势一致，拟合值与观测值重合程度高，模型的MAPE值为59.23%，观测值和拟合值都在预测值的95%置信区间内，绝对误差和相对误差均较小，说明该模型的预测精度较高，可以对未来进行很好地跟踪和预测。用该模型对青州市2018年各月布鲁氏菌病的新发病人数进行预测，并与实际月发病数进行比较，结果见表3。表3显示，实际值与预测值的平均误差为1.1，预测趋势与实际趋势基本一致，模型的预测效果较好。

图5 2011—2017年青州市布鲁氏菌病实际值与模型预测值比较

表3 2018年青州市布鲁氏菌病实际月发病数与模型预测月发病数比较

3 讨论

布鲁氏菌病在青州市的流行可追溯到1965年[9-10]，对青州市人群的健康影响较大。ARIMA模型是一种著名的时间序列统计模型,不需要分析影响疾病的各种因素，将各种影响因素的综合效应统一蕴涵于时间变量之中，研究时间变量的历史数据，挖掘出其中的规律以准确预测未来的发展变化趋势的数理统计模型[11]。

本文采用2011—2017年青州市布鲁氏菌病发病资料经过序列平稳化、模型识别和参数估计等一系列过程构建了ARIMA(0,1,1)(1,1,0)12预测模型，得到2018年布鲁氏菌病各月发病数预测值，并与实际值比较。预测结果显示布鲁氏菌病预测值和实际值的动态趋势基本一致，实际值均落入预测值的95%置信区间内,表明ARIMA(0,1,1)(1,1,0)12模型拟合效果较好，精度较高，体现出模型良好的实用性和应用价值[12]。因此ARIMA模型可以为传染病发病率进行早期预警，为传染病防控工作提供科学依据，具体可根据传染病既往的发展变化规律判断暴发或者流行的可能性，如果实际值在预测值95%置信区间范围内波动，表明疫情基本正常;如果超出预测值95%置信区间, 提示传染病暴发或流行的可能性较大[13]。

尽管ARIMA模型兼具回归分析和移动平均的优点，适用范围较因果回归分析法等分析方法广，研究过程相对简化、经济、适用，但是ARIMA模型是假定时间序列未来发展模式是由过去的惯性趋势发展来的，因此只适用于短期预测。而且ARIMA模型对数据的要求较高，要求时间序列在平稳的前提下要有30个以上数据。另外当实际情况复杂，特别是在采取了干预措施(如加强管理)或外部环境等因素发生较大改变时，模型的建立相对比较困难，应根据实际情况综合考虑预测结果[14]。可以在所建模型中不断增加样本量对模型进行修正甚至重新拟合，进一步提高模型的预测精度。

综上所述，本研究建立的ARIMA(0,1,1)(1,1,0)12模型能反映青州市布鲁氏菌病时间分布特征和发展趋势，预测效果较好，可以为疫情监测和疾病防控工作提供参考依据。