2016, Vol. 37
2. 第二军医大学基础医学部数理教研室, 上海 200433;
3. 解放军309医院全军结核病研究所, 北京 100091
2. Department of Mathematics & Physics, College of Basic Medical Sciences, Second Military Medical University, Shanghai 200433, China;
3. Institute for Tuberculosis Research, No. 309 Hospital of PLA, Beijing 100091, China
世界卫生组织(WHO)于2015年10月28日发布的《2015年全球结核病报告》中指出,2014年结核病在全球范围夺去了150万人的生命，仍然是最严重的公共卫生威胁之一[1]；WHO估算我国2014年新发肺结核人数为93万，仅次于印度和印度尼西亚位，居全球第3位[1]。虽然我国近年来逐步形成了由点到面、由局部到整体、由基层到中央的网络监测系统，但由于我国地域广阔、监管机构复杂、监测资料不完善等因素，缺少有效的预测预警分析手段，未能有效掌控肺结核病的发病趋势，使得该病防治工作处于被动的局面[2-3]。建立最优肺结核发病预测模型是正确预测肺结核发病水平、合理分配卫生资源和持续有效地开展肺结核病预警工作的重要前提。近年一些学者利用传染病监测数据建立自回归移动平均(autoregressive integrated moving average,ARIMA)预测疫情的发病率，取得了理想的预测效果[4-5]。本研究在时间序列分析基础上，采用乘积季节ARIMA(seasonal ARIMA,SARIMA)模型对我国2004年1月至2013年12月的肺结核逐月发病率进行拟合和预测，以期证明逼近准确的肺结核发病率预测可以量化对该病防治的决策。
1 资料和方法 1.1 资料来源原始数据来源于中华人民共和国国家卫生和计划生育委员会发布的2004年1月至2013年12月的全国法定传染病疫情概况以及国家统计局发布的2004 年至2013年人口统计资料。
1.2 研究内容通过研究发现我国肺结核月发病率序列的自相关系数明显有一个周期步长为12个月的强相关性，因此尝试采用时间序列分析法，对我国2004年1月至2012年12月的肺结核发病率原序列做周期步长为12个月的一阶季节差分变换，建立多个乘积SARIMA(p,d,q)×(P,D,Q)s模型，对肺结核逐月发病率进行拟合，通过时间序列特征分析、模型识别、参数估计及检验、模型诊断筛选出最优预测模型，并选取2013年1月至12月肺结核发病率数据以评价模型的预测性能。
1.3 研究方法 1.3.1 理论与模型介绍如果一个时间序列{Yt=yt-μ，t=1,2,3,…}满足如下模型[6]：
|
(1) |
则称其为季节周期s的乘积SARIMA(p,d,q)×(P,D,Q)s模型。其中：
|
(2) |
分别是p阶自回归过程AR(p)的自回归算子、q阶移动平均过程MA(q)的移动平均算子、P阶季节自回归算子、Q阶季节移动平均算子。
1.3.2 序列的平稳性检验、模型的识别与参数估计只有对于平稳时间序列，才可以使用自回归过程和移动平均过程。对于非平稳时间序列，可以尝试通过差分变换将其变为平稳时间序列。表现在上述模型中，就是差分的阶数d、季节周期s和季节差分的阶数D的选择问题。序列的平稳性检验通常有以下6种：ADF检验、DFGLS检验、PP检验、KPSS检验、ERS检验、NP检验。至于自回归的阶数p、季节自回归的阶数P、移动平均的阶数q以及季节移动平均的阶数Q的选择，一般是通过观察样本的相关图和偏相关图，根据其特点，初步为模型定阶。模型的定阶可能需要反复尝试并根据模型检验和拟合的结果进行选择，并没有一个普遍适用的做法[7]。初步确定模型的阶数后，使用EViews 7.0.0.1软件对模型参数进行估计。
1.3.3 模型的检验与诊断检验的内容主要包括：一是模型参数估计时对参数的检验；二是对模型的残差序列进行白噪声检验，以确保残差序列中不再包含还可以改进模型估计的有用信息。
1.3.4 模型的预测经过检验的模型，可以用于预测。对于SARIMA模型，预测一般分为样本内预测和样本外预测[8]。样本内预测又可以分为静态预测和动态预测。静态预测是指在预测时如果需要用到滞后期的数据，用真实数据代入；动态预测是指在预测时如果需要用到滞后期的数据，则用先前的预测数据代入，这相当于只用序列一开始的几个数据，预测其后所有的结果。因此，与静态预测相比，动态预测的效果通常要差。样本外预测，如果是多期预测，则只能用动态预测，如果是用静态预测，最多只能预测样本外的第一个时期。对模型预测结果的评估主要依据以下4个指标：误差均方根(RMSE)、平均绝对误差(MAE)、平均相对误差绝对值(MAPE)、Theil不等系数(TIC)。
本研究将全部观测分为两部分，一部分作为样本内数据(2004年1月至2012年12月的肺结核发病率)，用于模型的参数估计、识别及诊断，另一部分用于样本外预测(2013年1月至12月肺结核发病率)，并与实际观测比较，以考察模型的实用性及预测性能。
1.4 统计学处理应用Excel记录数据资料和模型预测结果；应用EViews 7.0.0.1进行SARIMA模型的参数估计、模型拟合及检验。
2 结 果 2.1 时间序列特征分析从样本内月发病率时间序列{IRt,t=1,2,…108}的自相关系数(图 1)明显可以看出该序列存在周期步长为12个月的强相关性，所以本研究尝试对原序列做周期步长为12个月的一阶季节差分变换，使样本内数据发病率序列D(IR,0,12)趋于平稳(图 2)。季节差分后序列平稳性检验结果显示，在平稳性1%水平下不能拒绝D(IR,0,12)是平稳序列。
|
图 1 样本内数据发病率序列的自相关系数 Fig 1 AC function of incidence rate sequence in the sample data |
|
图 2 原序列经季节差分后的自相关系数 Fig 2 AC function of original sequence with seasonal difference |
2.2 模型识别
对平稳化处理后的发病率序列建立乘积SARIMA(p,d,q)×(P,D,Q)s模型，根据样本相关图和偏相关图的特点反复试验，为模型定阶。
(1) 做周期步长为12个月的一阶季节差分变换，观察季节差分后的自相关系数图中AC值和PAC值(图 2)，发现图形的时点基本都在可信区间内，因此根据季节差分后自相关系数图的特点，初步确定模型形式为SARIMA(p,0,q)×(P,1,Q)12和SARIMA(p,1,q)×(P,1,Q)12。
(2) 一般情况下，p、q以及P、Q的识别较难，超过2阶的情况很少[9]，因此，通过观察季节差分后的自相关系数图中AC和PAC，初步确定P=2，q=1[10] 。p、q及P、Q采取0、1、2 组成的不同参数组合从低阶到高阶对模型进行反复调试和检验，确定3个备选模型：SARIMA(2,0,2)×(0,1,1)12；SARIMA(1,0,1)×(0,1,1)12；SARIMA(0,1,1)×(0,1,1)12。最后根据模型参数估计、残差序列的白噪声检验结果及拟合优度进行综合判断，筛选出最优预测模型。
2.3 模型参数估计及检验(1) 用EViews 7.0.0.1分析样本内数据后实现备选模型的参数估计(表 1)，由此可见备选模型的参数估计均具有统计学意义。(2) 参数估计之后，对模型的残差序列进行白噪声检验，以确保残差序列中不再包含还可以改进模型估计的有用信息(表 2)。根据各滞后期Q-Stat的P值，检验结果认为残差之间不存在相关性，即备选模型的残差序列是白噪声序列。(3) 备选模型拟合优度统计量比较显示，备选模型SARIMA(2,0,2)×(0,1,1)12的赤池信息准则(AIC)和Schwarz贝叶斯信息准则(SBC)值相对较小(表 3)，说明该模型与序列的拟合度最好。
|
表 1 备选模型参数估计及检验 Tab 1 Parameter estimation and test of alternative models |
|
|
表 2 备选模型残差序列的白噪声检验 Tab 2 White noise test for residual series of alternative models |
|
|
表 3 备选模型拟合优度统计量 Tab 3 Goodness of fit statistics of alternative models |
2.4 模型诊断
经过对备选模型的定阶、参数估计、残差序列白噪声检验及采用AIC和SBC衡量模型与序列的拟合度比较，最终选择我国肺结核发病率的预测模型为SARIMA(2,0,2)×(0,1,1)12，模型表达式为：
|
(3) |
其中IRt是发病率，▽12表示滞后期为12的一阶差分，▽12IRt=IRt-IRt-12，Ls是滞后算子，Lsεt=εt-s
2.5 模型验证与预测模型验证：由图 3可见，SARIMA(2,0,2)×(0,1,1)12模型样本内静态预测、样本内动态预测与肺结核实际发病数的拟合程度较好，实际值均在预测区间(95%CI)以内，且预测趋势与实际趋势基本吻合，符合本次建模要求。
|
图 3 SARIMA(2,0,2)×(0,1,1)12模型样本内静态预测拟合图(A)和动态预测拟合图(B) Fig 3 The fitting chart of static forecast (A) and dynamic forecast (B) in the sample by SARIMA(2,0,2)×(0,1,1)12 model |
模型预测：利用建立的SARIMA(2,0,2)×(0,1,1)12模型，预测2013年1月至12月发病率。由图 4可见，2013年1月至12月发病率真实值均在该模型预测区间(95%CI)以内，且预测趋势与实际趋势基本吻合。由表 4可见该模型预测值误差绝对值在0.036 673~1.239 393，误差绝对率在0.408 454%~16.179%。平均误差绝对值为0.416 992，平均误差绝对率为5.350 8%，提示该模型具有较好的预测性能。
|
图 4 SARIMA(2,0,2)×(0,1,1)12模型样本外动态预测拟合图 Fig 4 The fitting chart of dynamic forecast outside the sample by SARIMA(2,0,2)×(0,1,1)12model |
|
|
表 4 SARIMA(2,0,2)×(0,1,1)12模型对2013年肺结核发病率的预测结果 Tab 4 SARIMA(2,0,2)×(0,1,1)12 model in prediction of pulmonary tuberculosis incidence of 2013 |
3 讨 论
建立最优预测模型，提高预测准确性，对于肺结核病防治工作意义重大。在肺结核流行趋势的预测中，数学模型起着极其重要的作用。目前有多种数学模型被应用于传染病预测，如灰色预测模型、马尔科夫链预测模型[11]、回归模型[12]、神经网络模型(ANN)[13]、ARIMA模型等[14]。其中ARIMA模型是对有时间性变动的序列提取季节趋势建立模型的方法，其精确度较高，是传染病时间序列预测模型中最重要的手段[15]。因此，对于肺结核发病率的预测，采用乘积SARIMA模型是较为理想的。
本研究基于样本内数据，通过时间序列特征分析、模型识别、模型参数估计及检验、模型诊断建立最优预测模型——乘积SARIMA(2,0,2)×(0,1,1)12模型。预测结果显示2013年1月至12月发病率真实值均在该模型预测区间(95%CI)以内，且预测趋势与实际趋势基本吻合，平均误差绝对值为0.416 992，平均误差绝对率为5.350 8%，提示该模型对肺结核发病趋势进行了较准确的跟踪和预测，能较好地模拟和预测肺结核发病率在时间序列上的变动趋势，可以为肺结核的控制和管理提供量化参考依据。
乘积SARIMA模型是在样本内数据发病率序列趋于平稳的基础上较好地拟合了历史数据，是精度较高的短期预测方法。当用于长期预测时，应根据监测数据调整模型的各项参数。若序列太短，则可靠性较差。在时间序列预测模型中，对于模型的定阶，没有固定的统一的规则可以严格遵循，因此，模型的定阶过程存在一定的主观性，需要摸索和试错。如果环境因素或其他影响发病率的外界因素发生较大的变化或突变，则原来时间序列的内部规律会受到冲击和破坏，预测模型的外推(样本外的预测)可能不再有效。
乘积SARIMA模型是从数据上反映肺结核发病率的发展趋势，可以为肺结核的防治工作提供数据参考依据，建议有关部门在制定具体肺结核防控策略时还应考虑其他综合因素的影响。
| [1] | World Health Organization. Global Tuberculosis Report 2015[R/OL].[2016-07-02].http://www.who.int/tb/publications/global_report/gtbr15_executive_summary_zh.pdf?ua=1. |
| [2] | 杨召, 叶中辉, 尤爱国, 郭奕瑞, 张肖肖, 梁淑英, 等. 乘积季节ARIMA模型在结核病发病率预测中应用[J]. 中国公共卫生 , 2013, 29 :469–472. |
| [3] | 牛成虎, 梅光辉, 石敏, 高红韦. 我国肺结核发病率的发展动向及预测研究[J]. 现代生物医学进展 , 2009, 9 :561–564. |
| [4] | HU W, TONG S, MENGERSEN K, CONNELL D. Weather variability and the incidence of cryptosporidiosis: comparison of time series poisson regression and SARIMA models[J]. Ann Epidemiol , 2007, 17 :679–688. DOI:10.1016/j.annepidem.2007.03.020 |
| [5] | GOMEZ-ELIPE A, OTERO A, VAN HERP M, AGUIRRE-JAIME A. Forecasting malaria incidence based on monthly case reports and environmental factors in Karuzi, Burundi, 1997-2003[J]. Malar J , 2007, 6 :129. DOI:10.1186/1475-2875-6-129 |
| [6] | 赵国庆. 经济分析中的时间序列模型[M]. 天津: 南开大学出版社, 2012 : 17 -22. |
| [7] | CRYER J D,CHAN K S. 时间序列分析及应用[M].2版.潘红宇,译.北京:机械工业出版社,2011:40-58. |
| [8] | 易丹辉. 数据分析与EViews应用[M]. 北京: 中国人民大学出版社, 2008 : 141 -148. |
| [9] | 温亮, 徐德忠, 林明和, 夏结来, 张治英, 苏永强. 应用时间序列模型预测疟区疟疾发病率[J]. 第四军医大学学报 , 2004, 25 :507–510. |
| [10] | 宇传华. SPSS与统计分析[M]. 北京: 电子工业出版社, 2007 : 577 -614. |
| [11] | UYS P W, VAN HELDEN P D, HARGROVE J W. Tuberculosis reinfection rate as a proportion of total infection rate correlates with the logarithm of the incidence rate: a mathematical model[J]. J R Soc Interface , 2009, 6 :11–15. DOI:10.1098/rsif.2008.0184 |
| [12] | SHILOVA M V, GLUMNAIA T V. [Prediction of the rate of tuberculosis mortality (calculation methodology)][J]. Probl Tuberk Bolezn Legk , 2006 (1) :22–28. |
| [13] | HAMDY K E, YOSRY A A, FARAG I Y. Prediction of hourly and daily diffuse fraction using neural network, as compared to linear regression models[J]. Energy , 2007, 32 :1513–1523. DOI:10.1016/j.energy.2006.10.010 |
| [14] |
胡晓媛, 吴娟, 孙庆文, 沙琨, 王玲玲, 李敏. ARIMA模型与GRNN模型对肺结核发病率预测的对比研究[J]. 第二军医大学学报 , 2016, 37 :115–119.
HUI X Y, WU J, SUN Q W, SHA K, WANG L L, LI M. Comparative study on ARIMA model and GRNN model for predicting the incidence of tuberculosis[J]. Acad J Sec Mil Med Univ , 2016, 37 :115–119. |
| [15] | 孟蕾, 王玉明. ARIMA模型在肺结核发病预测中的应用[J]. 中国卫生统计 , 2010, 27 :507–509. |