近年来，手足口病的暴发或流行，引起儿童死亡，受到越来越多的关注^[1-2]。2008年5月将手足口病纳入丙类法定传染病管理以来，陕西省报告的发病数和死亡数逐年上升。针对陕西省手足口病的研究主要集中在病原学及流行特征分析^[3-4]，本研究以R软件为基础，应用自回归求和移动平均（autoregressive integrated moving average，ARIMA）乘积季节模型，对陕西省手足口病发病人数进行分析预测，以了解陕西省手足口病发病人数的变化规律和发展趋势。

资料与方法

1. 数据来源：2009－2015年陕西省各级各类医疗卫生机构通过《传染病报告管理信息系统》报告的手足口病发病人数。2009年1月至2015年6月手足口病发病人数数据用于模型拟合，2015年7月至2015年12月手足口病发病人数数据用于模型预测效果评价。

2. 研究方法：ARIMA模型是一种基于时间序列分析和预测的方法。随机季节模型与ARIMA模型结合即为ARIMA季节模型，一般表示为ARIMA（p，d，q）×（P，D，Q）s，各参数的意义：非季节自回归阶数p、非季节差分阶数d、非季节移动平均阶数q，季节自回归阶数P、季节差分阶数D和季节移动平均阶数Q，季节长度s。当P＝D＝Q＝0时，该模型为ARIMA模型。ARIMA乘积季节模型能综合考虑季节、趋势和随机干扰等因素，对时间序列预测效果较好^[5]。此次研究采用R-3.2.2软件进行建模及数据处理分析。建模过程：①序列平稳：通过“tseries”软件包中的“ts”函数将2009年1月至2015年6月手足口发病人数定义为以月为单位的时间序列，利用“decompose”函数绘制并分解时间序列，观察时间序列的趋势及周期变化，然后通过“diff”函数对序列进行差分，使之成为零均值、无明显趋势变化的平稳序列。②模型识别：对平稳时间序列，分别通过“acf”和“pacf”函数观察序列的自相关系数（autocorrelation function，ACF）和偏自相关系数（partial autocorrelation function，PACF）初步确定p，d，q和P，D，Q的阶值。若自相关函数为拖尾，且偏相关函数为截尾，适合AR模型；若自相关函数为截尾，且偏相关函数为拖尾，适合MA模型；若自相关函数和偏相关函数均为拖尾，则序列适合ARIMA模型。参数P、Q超过2阶的情况不多见，可分别取0、1、2并由低阶到高阶逐个试验，根据模型的拟合优度等指标进行判断，以确定这两个参数。s根据疾病的背景知识获得。③参数估计：ARIMA乘积季节模型运用最大似然估计（maximum likelihood estimation，MLE）或最小二乘估计（least squares estimation，LSE）法，根据p，d，q和P，D，Q的阶值计算出自回归和移动平均系数。④模型诊断：通过“Box.test”函数计算Box-Ljung统计量，如果统计量差异无统计学意义，可认为残差为白噪声，若残差序列不是白噪声，则需对模型进行改进。若残差序列满足白噪声过程，通过AIC等指标挑选最优模型，AIC越小，模型的拟合程度越好。此外，通过比较2015年7－12月手足口病实际发病人数和模型拟合发病人数，计算各月绝对误差和相对误差水平来评价模型的拟合效果。⑤模型预测：以上步骤建立2009年1月至2015年6月陕西省手足口病发病人数模型，通过“forecast.Arima”函数模型预测2016－2017年发病人数。

3. 统计学分析：本研究采用Excel软件建立数据库，用R-3.2.2中“tseries”和“forecast”软件包进行数据处理和建立模型。R软件是由奥克兰大学的Robert Gentlman和Ross Ihaka及其他志愿人员开发的统计分析软件，可用来探索数据、统计分析及做图等。

结 果

1. 发病人数变化趋势：2009年1月至2015年6月陕西省手足口病年平均发病人数为52 631.14例，月平均发病人数4 385.93例，其中发病人数最多的是2015年5月，共发病18 074例。由时间序列可看出陕西省7年间手足口病发病人数存在一定的周期性，将时间序列用“decompose”语句按总体趋势、季节趋势及随机误差分解后可看出，总体趋势除在2011年出现下降外，2009－2015年手足口病发病人数呈整体上升趋势，提示该序列为非平稳时间序列。此外，手足口病存在明显的季节性，春夏之交出现高峰，两头下降，随机误差保持在一定水平（图 1）。

2. ARIMA乘积季节模型构建：

（1）序列平稳化：由图 1可知时间序列呈现明显的季节变化。用“diff”语句进行一阶12步季节差分消除时间序列趋势和季节影响，得到了一个近似平稳的随机序列（图 2），符合ARIMA乘积季节模型对平稳性的要求。

（2）模型识别：根据差分变换的次数，初步确定以12个月为周期的ARIMA（p，1，q）（P，1，Q）12乘积季节模型，分别用“acf”和“pacf”语句画出自相关图（图 3）和偏自相关图（图 4）。由图可知，自相关系数在0阶截尾，12阶、24阶有明显波动，说明差分后序列中含有显著地季节效应，偏自相关系数在2阶截尾，同时在滞后12阶显著大于可信区间范围，也显示出明显的季节周期，考虑模型的AR因子为（1－θ₁B－θ₂B²）ε_t。根据上述特征判断p＝2，q＝0，但P和Q的判断比较困难，根据《现代医学统计学》（方积亁，陆盈. 人民卫生出版社，2002），P和Q超过2阶的情况很少见，可分别取0、1、2逐个试验，根据参数及模型总体的显著性、拟合优度等指标进行比较，选择最佳模型。

（3）参数估计和模型诊断：拟合的3个备选模型参数估计值及检验结果见表 1。残差白噪声检验和拟合优度检验用来评估各个备选模型对时间序列内在综合信息的提取程度及各个备选模型之间的优劣比较。用“Box.test”语句“Ljung-Box”统计量对3个备选模型的残差值进行检验，残差均为白噪声，由AIC统计量可知，模型1拟合效果最好（表 1）。综合以上结论，选择模型ARIMA（2，1，0）×（1，1，0）12是合适的。

3. 模型评价：利用ARIMA（2，1，0）×（1，1，0）12乘积季节模型对陕西省2015年7－12月手足口病发病人数进行预测（表 2）。该模型预测值的绝对误差、相对误差范围分别为146.520～1 843.137和0.053～0.188，绝对误差的平均531.535，相对误差平均0.114，提示ARIMA（2，0，0）×（1，1，0）₁₂模型具有较佳的预测性能。

4. 模型预测：用“forecast”函数包中的“plot.forecast”过程，绘制模型拟合的陕西省2009年1月至2015年12月手足口病发病人数（虚线）及2016－2017年手足口病发病人数（实线）序列图（图 5）。提示2016年和2017年陕西省手足口病发病水平与2015年接近。

讨 论

自2009年开展监测以来，手足口病发病人数不断呈上升趋势，而发病人数的预测对该病的防治工作具有重要指导意义。数学模型预测是进行手足口病预警及控制的较好手段。手足口病预测模型主要建立在该病疫情发展的基础上。目前手足口病的预测已有多种数学模式^[6-8]，但均需要传播环节或敏感人群的分布参数。一旦参数发生变化或无法获得时，也将无法做出准确预测。此外，实际工作中，仍有很多影响手足口病的因素未被考虑到模型中，从而影响到结果的准确性和可延性。而时间序列预测不需考虑实际参数，而是在长期的实际情况下达到平衡状态时做出预测^[9-10]，手足口病的发病人数序列具有非线性和季节性，给P，Q的取值带来困难。

本研究用开源的R软件以及相关软件包^[11]，通过ARIMA乘积季节模型分析2009年1月至2015年6月陕西省手足口病的时间序列数据，探讨该病的发病规律，并对发病人数作出预测。该模型可为陕西省手足口病预防措施的制定提供一定依据。作为一种数据处理方法，时间序列预测方法与其他常用的预测方法一样，主要从数理上反映疾病的统计规律，作为疾病防治的决策还须考虑其他因素对预测结果的影响。另外，本文提出的ARIMA乘积季节模型能预测不同时间的手足口病发病人数，但对短期预测还是长期预测较为准确，还需进一步研究和后续数据的支持。