张顺先 邱磊 张少言 李翠 胡骏 田黎明 鹿振辉
结核病是由结核分枝杆菌引起的全球最流行的传染性疾病之一,可发生于全身多脏器,以肺部受累最为常见[1],是世界上因单一病原体感染造成死亡最多的传染病[2-3]。2019年全球新发肺结核患者超过1000万例,肺结核相关的死亡患者超过120万例[2]。近年来中国作为全球30个结核病高负担国家之一[4],肺结核疫情具有患病率下降缓慢、疾病负担巨大和地区不平衡等特点,防治任务艰巨[1, 5],准确掌握肺结核新发患者例数及流行特征,是开展肺结核相关健康教育、检测技能培训和配置相关卫生资源的前提和基础。伴随计算机软件的进步,很多数学模型因其能够较好地揭示传染病随时间发展的动态变化和预期结果相对准确等特点被广泛应用于传染病疫情的预测和研究,为疾病的防控提供积极参考[6-7]。
自回归移动平均(autoregressive integrated moving average,ARIMA)模型从序列自相关的角度揭示了传染病的时间序列发展规律,既能解决数据间的自相关问题,又能充分利用其中的季节效应特征拟合疾病发生的实际情况,短期预测效果较好[8-9],被广泛应用于传染病的预测预警。本研究采用该模型对我国(不含我国香港、台湾、澳门地区,下同)肺结核月报告患者例数进行时间序列发展规律研究,以预测中国肺结核月报告患者例数,为结核病相关药品的生产、采购和就诊入院等卫生服务的供给预算提供科学参考。
一、资料来源
通过中国疾病预防控制中心主办的《疾病监测》杂志公布的中国每月甲、乙、丙类传染病疫情动态简介,搜集2006年1月至2019年8月全国肺结核月报告患者例数(不含中国港澳台地区)。采用SPSS 26.0软件对2006年1月至2018年12月的全国肺结核月报告患者例数进行拟合ARIMA备用模型,再以2019年1—8月的数据评价备用模型的预测效果和筛选最优模型。疫情动态数据由中国各级各类医疗卫生机构通过《传染病报告管理信息系统》在线报告、后期复核而产生。
二、ARIMA建模
1.时间序列的建立、特征分析及平稳化处理:建立时间序列,然后通过自相关系数函数图(autocorrelation function,ACF)和偏相关系数函数图(partial autocorrelation function,PACF)分析序列的平稳性,如果序列不平稳,将其通过自然对数转换、一般差分和(或)季节差分等方法处理数据[10]。
2.模型识别:根据ACF图、PACF图和SPSS 26.0统计学软件初步拟合的结果,尝试对模型进行初步识别和定阶。预先考虑使用ARIMA(p,d,q)或ARIMA(p,d,q)×(P,D,Q)模型,结合以往经验,p和q一般不会超过3阶,d不会超过2阶,P、D和Q不会超过2阶,p、q、P和Q可依次取0、1和2并由低阶到高阶逐个尝试,同时d和D依次取0和 1进行尝试[10-11]。
3.备选模型的整体诊断和参数判读:通过非线性最小二乘法估计模型参数,然后对模型残差进行白噪声(Ljung-BoxQ)检验,以判断模型整体的合理性。然后再根据模型整体的检验指标(选择P>0.05的统计量,且Ljung-BoxQ值越小越好)、整体模型是否简洁[即p、d、q和(或)P、D、Q的阶越低越简洁]、ARIMA模型各参数[自回归法(autoregressive,AR),平均移动法(moving average,MA),季节自回归法(seasonal autoregressive,SAR),季节移动平均法(seasonal moving average,SMA)]是否有意义、整体模型的平稳决定系数(stationary R-square,R2;越大越平稳)、整体模型的标准化贝叶斯信息准则值(normalized bayesian information criterion,NBIC;越小越好),及整体模型的均方根误差(root mean square error,RMSE;越小越好)来筛选备选模型[7, 12]。
图1 2006—2018年全国肺结核月报告患者例数时间序列图
4.最优模型的确定:在满足模型成立的基本条件、完成备选模型的选择后,从预测值相对误差越小模型越优的原则出发[13],应用不同的备选模型去预测2019年1—8月的全国肺结核月报告患者例数,再与2019年1—8月实际报告的患者例数进行比对,选择预测平均相对误差最小的模型为最优模型。
5.模型的预测应用:使用最优模型预测2019年9月至2020年12月我国每月的肺结核新发患者例数。
三、 统计学处理
采用Excel 2019软件对2006年1月至2019年8月的全国肺结核月报告患者例数进行整理,采用SPSS 26.0软件进行ARIMA模型预测分析,通过比较预测数据和实际数据的相对误差来确定最优模型,并进行2019年9月至2020年12月数据的预测应用。
一、中国肺结核月报告患者例数的发病趋势分析
2006年1月至2018年12月中国肺结核报告患者总计为16 905 732例,平均年报告患者例数为1 300 441 例,从2006—2008年,报告患者例数缓慢上升(2006年为1 454 232例、2007年为1 389 072例、2008年为1 532 472例),从2009年开始年报告患者例数逐步下降。2006年1月至2018年12月中国肺结核平均月报告患者例数为108 370例,其中,2015年2月报告患者例数最少(75 541例),2008年3月报告患者例数最多(156 679例)。中国肺结核报告患者例数总体上呈现先短暂上升再缓慢下降的流行特征(图1),即非平稳序列,同时也呈现出一定的季节特征,每年2—4月为流行高峰、11月到次年1月为流行低谷。
二、模型的参数估计和模型诊断
以2006年1月至2018年12月中国肺结核月报告患者例数建立时间序列,发现其ACF图和PACF图的相关系数在延迟数目12[即“滞后数(lags)”]时尚未落入95%的置信区间内,说明序列(连续12个月的数据)不稳定,原始数据可能在12个月的周期上有内在关联,遂将数据进行自然对数转化和1阶差分处理,需建立模型为ARIMA(p,d,q)或ARIMA(p,d,q)×(P,D,Q)。在模型总体的Ljung-BoxQ值所对应的P值均>0.05的情况下,再根据模型是否简洁、模型各参数差异是否均有统计学意义(要求各参数的t值所对应的P值均<0.05)筛选出12个基本模型(表1)。然后再从这些基本模型中筛选出R2最大的模型ARIMA(1,0,1)(0,1,1)12,R2=0.707],RMSE最小的模型
表1 根据2006—2018年中国肺结核月报告患者例数建立的ARIMA模型
注AR:自回归法;MA:平均移动法; SAR:季节自回归法;SMA:季节移动平均法;R2:整体模型的平稳决定系数;NBIC:整体模型的标准化贝叶斯信息准则值;RMSE:整体模型的均方根误差
ARIMA(0,1,2)(0,1,1)12,RMSE=9147.85]、NBIC最小的模型[ARIMA(0,1,1)(0,1,1)12,NBIC=18.355]、Ljung-BoxQ值最小的模型[ARIMA(1,1,1)(0,1,1)12,Ljung-BoxQ=8.797]作为备用模型。
三、最优模型的确定
从预测值相对误差越小模型越优的原则出发确定最优模型。分别以备用模型ARIMA(0,1,1)(0,1,1)12、ARIMA(0,1,2)(0,1,1)12、ARIMA(0,1,1)(0,1,1)12和ARIMA(1,1,1)(0,1,1)12预测2019年1—8月中国肺结核月报告患者例数(表2),并与实际的月报告患者例数进行比较,发现ARIMA(0,1,1)(0,1,1)12模型的预测平均相对误差最小(0.55%),预测效果较好;其余3种模型的预测平均相对误差均较大(分别为3.01%、1.16%和0.58%),预测效果相对不理想。另外,从预测值最大相对误差的角度讲,ARIMA(1,0,1)(1,0,1)12的最大相对误差为10.88%、ARIMA(1,0,1)(1,0,1)12的最大相对误差为8.56%、ARIMA(1,1,2)(0,1,1)12的最大相对误差为9.58%,均高于ARIMA(0,1,1)(0,1,1)12的最大相对误差(8.17%),故最优模型确定为ARIMA(0,1,1)(0,1,1)12。
四、模型的预测应用
通过ARIMA(0,1,1)(0,1,1)12模型预测2019年9月至2020年12月中国肺结核月报告患者例数(图2),结果显示2019年9—12月预测患者例数分别为84 399(95%CI:98 844,71 586)例、79 928(95%CI:93 708,67 716)例、82 551(95%CI:96 884,69 858)例、82 598(95%CI:97 042,69 818)例;2020年1—12月预测患者例数分别为78 381(95%CI:92 488,65 941)例、75 614(95%CI:89 338,63 524)例、101 847(95%CI:120 486,85 444)例、95 075(95%CI:112 616,79 653)例、94 289(95%CI:111 825,78 887)例、89 862(95%CI:106 707,75 081)例、88 964(95%CI:105 770,74 231)例、84 980(95%CI:101 157,70 811)例、81 188(95%CI:96 761,67 561)例、76 873(95%CI:91 730,63 886)例、79 380(95%CI:94 834,65 882)例、79 410(95%CI:94 984,65 821)例(图2)。预计2020年新发肺结核患者约1 025 863例,平均每月85 489例。
近来年,ARIMA模型、灰色模型和神经网络等数理统计学因其能相对准确的探索传染病的发生发展规律,在传染病预测研究领域的应用非常活跃,为传染病的有效防控提供了科学参考和重要依据。ARIMA模型之所以能在疾病预测领域中广泛应用,因其既吸收了传统回归分析的优点又发挥了移动平均的长处,可以将诸多影响传染病发生发展的复杂因素(如人口流动、气温和湿度等因素)综合统一并蕴含于时间变量中,借助趋势变化、周期变化、随机干扰和模型参数进行量化表达,是一种适用范围广、实用性强、预测精确度较高、预测效果优于回归模型的预测方法。但从传染病的长期趋势来看,疫苗的出现和推广会极大地改变相关传染病的流行特征,如果采用ARIMA进行长期预测会出现一定的偏差,故ARIMA只适合短期预测,不适合长期预测[14]。另外,由于用ARIMA进行预测时一般需要10个以上的周期性变化的数据,故在急性传染病的早期阶段、未出现疾病流行的周期性变化趋势之前,ARIMA不能做预测,如2019年年底开始在全球流行的新型冠状病毒肺炎疾病(COVID-19),直至2020年3月30日仍不具备一个完整的周期性规律变化,如果以这段时间每天的报告患者例数进行ARIMA 预测,结果将出现很大的偏差。而对那些已出现过多个周期性趋势的慢性传染病(如肺结核),因其致病因素、影响因素、防控措施和策略(如疫苗等)不会在短时间内出现大的改变,导致ARIMA 对这类疾病的短期预测效果较好。
表2 2019年1—8月中国肺结核患者报告例数预测值与实际值比较情况
注平均相对误差为2019年1—8月相对误差之和/8×100%
图2 最优模型ARIMA(0,1,1)(0,1,1)12预测2019年9月至2020年12月中国肺结核月报告患者例数拟合图
肺结核是我国常见的一种对人群威胁较大的呼吸道传染病。本研究根据2006年1月至2018年12月我国传染病监测系统中肺结核报告病例数的变化趋势,发现中国肺结核新发患者例数总体上呈缓慢下降的趋势,可能与这些年中国肺结核大力推行短程督导化疗方案有关,同时也与大力推广肺结核患者的早发现、早诊断和早治疗诊治措施有关[15-16]。由于高速铁路网的迅速发展加剧了人口迁移和流动,以及结核分枝杆菌耐药性的出现和逐渐加重使得肺结核新发患者例数下降缓慢[17-18]。本研究发现中国肺结核发病总体在时间分布上具有鲜明的季节特征和周期变化趋势,我国每年肺结核的发病高峰为2—4月,这可能与春节期间大量的结核病患者不按规定进行自我隔离治疗,不采取戴口罩等防护措施,频繁的流动和参加聚集性活动,造成较大范围的传播有关[17-18]。同时,2—4月是我国的冬春季,天气寒冷,空气污染指数相对较高,不利于家庭开窗通风,增加了呼吸道感染(包括肺结核)的发生和发展[3, 19]。因此,提高肺结核患者预防传染病的自我意识,使其主动在家隔离、治疗并佩戴口罩外出;同时,建议家庭每天开窗通风,可在一定程度上降低发生肺结核的风险。
本研究应用ARIMA模型对我国近年来肺结核新发患者例数进行分析,在掌握肺结核发病规律和流行特征的同时,探索预测效果较好的模型。准确预测肺结核的发病趋势,不仅可评估防控措施的实施效果,同时对生产和采购抗结核相关药品、诊断试剂、耗材和卫生服务、合理分配卫生资源、制定最优的控制策略和措施具有重要意义,有助于提升肺结核防控的预见性和主动性。本研究通过预测发现2020年我国新发肺结核患者约1 025 863例,平均每月85 489例,这和近几年中国肺结核流行强度基本一致[12],但稍高于WHO的预测结果(90万新发患者)[1],这可能与ARIMA更适合短期预测有关,ARIMA进行长期预测可能会高估结果。因此,我们认为在没有其他因素干扰的情况下,2020年结核病防控相关的医疗资源、防控措施和力度仍需保持近2年的强度,不能降低各种相关资源的投入。
综上所述,ARIMA模型对中国肺结核新发患者例数预测效果较好,可为肺结核的防治提供科学参考。但模型的建立和预测应用是个动态过程,需不断根据积累的数据进行调整,从而提高预测精度。