手足口病是以肠道病毒enterovirus 71(EV71)和Coxsackievirus A 16(CoxA16)为主的多种肠道病毒引起的急性传染性疾病，临床表现以手、足、口腔等部位出现斑丘疹、疱疹为主要特征，多发于＜3岁学龄前儿童^[1]。手足口病于2008年5月被纳入国家丙类法定报告传染病，2009—2013年中国年均报告发病数为170多万，且有逐年上升趋势^[2]。据报道，湖南省手足口病疫情一直呈高发态势,2010和2012年，不仅发病数占全国前5位，重症病例、死亡病例报告数也是全国重点地区^[3]，其预防和控制问题事关重大。传染病研究中的数学建模分析结果能提供一定的理论基础和概念，通过数学模型可发现传染病的传播机制，预测传染病的流行趋势^[4]。季节性差分自回归移动平均模型(seasonal autoregressive integrated moving average model，SARIMA)属于时间序列模型的一种，是将非平稳时间序列平稳化，然后将因变量仅对它的滞后值以及随机误差项的现值和滞后值进行回归所建立的模型，用于对过去的变化趋势延续到未来的数据资料进行分析^[5]，是目前国内常用的传染病预测模型之一，在手足口病预测领域也有广泛应用^[6]。本研究拟采用SARIMA模型对湖南省手足口病疫情进行预测，为湖南省手足口病的防控提供科学依据。

资料来源于“中国疾病预防控制信息系统”中2008年5月—2014年7月按发病日期统计的湖南省各月手足口病报告月发病率，包括临床诊断病例和实验室诊断病例。

采用SPSS 18.0软件进行统计分析。以2008年5月—2013年12月的手足口病月发病率数据为基础，在“Define Date”模块和“ForeCasting”模块中分析，通过创建时间序列、数据平稳化处理、模型的识别、模型的参数估计和白噪声检验等步骤建立模型^{[7, 8]}，对2014年1月—2015年7月的发病率进行预测，并以2014年1-7月的实际数据进行验证。SARIMA模型建立在差分自回归移动平均模型(autoregressive integrated moving average model，ARIMA)^[9]基础上，增加了对季节性和周期性的分析。SARIMA模型的数学表达式为：φ(B)(1-B)^dφ(B^S)(1-B^S)DYt=ψ(B)Φ(B^S)εt。其中，Y_t=φ₁Y_t-1+φ₂Y_t-2+…+φ_pY_t-p+ε_t-ψ₁ε_t-1-ψ₂ε_t-2-…-ψ_qε_t-q；φ(B)=1-φ₁B¹-φ₂B²-…-φ_pB^p，φ(B^S)=1-φ₁B^S,1-φ₂B^S,2-…-φ_pB^S,P；ψ(B)=1-ψ1B-ψ₂B²-…-ψ_pB^q,Φ(B^S)=1-Φ₁B^S,1-Φ₂B^S,2-…-Φ_QB^S,Q。B是后移算子，p是AR(p)项的阶数，q是MA(q)项的阶数,d是简单差分的阶数；ε_t是时间序列在t期的误差或偏差，ε_t-p是时间序列在t期的误差或偏差；P是SAR(P)项的阶数，Q是SMA(Q)项的阶数，D是季节性差分的阶数。SARIMA模型可以表示为SARIMA(p,d,q)×(P,D,Q)_S，包含非季节性部分(p,d,q)和季节性部分(P,D,Q)_S，非季节性部分即ARIMA(p,d,q)^[10]。

以2008年5月为第1个个案时间定义日期，创建时间序列，生成2008年5月—2013年12月湖南省手足口病疫情月发病率的序列图。图 1显示，湖南省手足口病疫情月发病情况存在明显的周期性和季节性，约12个月为1个周期,每年5—7月出现第1次大高峰，11月出现第2次小高峰，属于不平稳序列。

图 1

Fig. 1

图 1 湖南省2008—2013年手足口病月发病率(1/10万)

为了降低方差以及消除季节性的影响，对原始时间序列进行1次平方根转换和1次季节性差分。图 2显示，进行1次季节性差分后的序列接近平稳，由此初步暂定SARIMA(p,d,q)×(P,D,Q)₁₂模型中的d和D分别为0和1。图 3显示，自相关系数(auto correlation function，ACF)呈缓慢余弦周期性衰减，拖尾，具有伪周期性，提示是平稳序列。图 4显示，偏自相关系数(partial auto correlation function，PACF)在滞后值lag=13处近似截尾。

图 2

Fig. 2

图 2 湖南省2008—2013年手足口病月发病率季节性差分序列图

图 3

Fig. 3

图 3 湖南省2008—2013年手足口病月发病率季节性差分ACF图

对比常见ARIMA模型自相关分析标准图，提示为AR过程，可以确定模型中非季节性部分中的p和q分别为1和0，模型可暂定为SARIMA(1,0,0)(P,1,Q)₁₂。对于季节性部分的P、Q，分别取值0、1和2进行逐个试验。表 1显示，9种备选模型的平稳R²、均方根误差和标准化的贝叶斯信息准则(Bayesian information criterion，BIC)值相差不大，其中SARIMA(1,0,0)(0,1,0)₁₂模型的Ljung-Box Q统计量的P<0.05，故可先排除。结合9种备选模型的参数估计表(表 2)可知，所有的参数均有统计学意义的模型有SARIMA(1,0,0)(0,1,0)₁₂(已排除)、SARIMA(1,0,0)(1,1,0)₁₂，故最终选择的模型为SARIMA(1,0,0)(1,1,0)₁₂，该模型的平稳的R²为0.779，均方根误差为0.653，标准化的BIC值为-0.637。

图 4

Fig. 4

图 4 湖南省2008—2013年手足口病月发病率季节性差分PACF图

表 1 (Table 1)

表 1 备选9种模型的模型拟合和统计 模型 模型拟合统计量 Ljung-Box Q 平稳的R2 R2值 均方根误差 标准化的BIC 统计量 df P值 (1,0,0)(0,1,0)12 0.627 0.728 0.840 -0.204 40.158 17 0.001 (1,0,0)(0,1,1)12 0.753 0.820 0.690 -0.528 13.168 16 0.660 (1,0,0)(0,1,2)12 0.758 0.824 0.689 -0.456 13.741 15 0.545 (1,0,0)(1,1,0)12 0.779 0.839 0.653 -0.637 15.420 16 0.494 (1,0,0)(1,1,1)12 0.785 0.844 0.649 -0.576 15.442 15 0.420 (1,0,0)(1,1,2)12 0.789 0.847 0.650 -0.502 12.789 14 0.543 (1,0,0)(2,1,0)12 0.782 0.842 0.654 -0.563 15.836 15 0.393 (1,0,0)(2,1,1)12 0.786 0.845 0.654 -0.491 14.522 14 0.412 (1,0,0)(2,1,2)12 0.786 0.845 0.661 -0.398 14.705 13 0.326

表 1 备选9种模型的模型拟合和统计

表 2 (Table 2)

表 2 9种备选模型的参数估计 模型 常数项 AR1 SAR1 SAR2 SMA1 SMA2 B值 t值 P值 B值 t值 P值 B值 t值 P值 B值 t值 P值 B值 t值 P值 B值 t值 P值 (1,0,0)(0,1,0)12 0.247 0.500 0.608 0.793 9.830 <0.001 (1,0,0)(0,1,1)12 0.487 2.959 0.005 0.770 8.922 <0.001 0.979 0.263 0.794 (1,0,0)(0,1,2)12 0.491 2.582 0.013 0.773 8.912 <0.001 0.760 4.458 <0.001 -0.179 -0.752 0.456 (1,0,0)(1,1,0)12 0.489 2.337 0.023 0.765 8.789 <0.001 -0.685 -6.229 <0.001 (1,0,0)(1,1,1)12 0.485 2.087 0.042 0.754 8.439 <0.001 -0.868 -5.477 <0.001 -0.362 -1.193 0.238 (1,0,0)(1,1,2)12 0.533 2.724 0.009 0.761 8.313 <0.001 -0.994 -4.356 0.001 -0.482 -0.321 0.749 0.394 0.457 0.650 (1,0,0)(2,1,0)12 0.477 2.044 0.046 0.761 8.654 <0.001 -0.585 -3.823 <0.001 0.156 0.795 0.431 (1,0,0)(2,1,1)12 0.511 2.281 0.027 0.759 8.341 <0.001 -1.386 -1.365 0.178 -0.392 -0.537 0.593 -0.907 -0.417 0.679 (1,0,0)(2,1,2)12 0.507 2.220 0.031 0.759 8.209 <0.001 -1.475 -2.800 0.007 -0.479 -1.125 0.266 -0.992 -0.604 0.549 -0.057 -0.138 0.891

表 2 9种备选模型的参数估计

对选定的模型进行白噪声检验，残差的ACF和PACF图见图 5。图 5显示，残差的ACF和PACF均在95%CI以内，且Ljung-Box Q=15.420，P=0.494，残差序列为白噪声，表明拟合的该模型合理。

图 5

Fig. 5

图 5 模型SARIMA(1,0,0,)(1,1,0)12残差序列ACF、PACF图

用选定的SARIMA(1,0,0)(1,1,0)₁₂模型对2014年1月-2015年8月的湖南省手足口病月发病率情况进行预测，其预测值和95% confidence interval(95%CI)见表 3、4。图 6显示，观测值与拟合值的重合度较高，基本均在拟合值的95%CI以内。预测2014年湖南省手足口病发病率第1次高峰在5月份，第2次高峰在11月份，2次高峰期的发病率均高于2013年同月份的发病率；2015年上半年高峰期也在5月份，高峰期的发病率低于2014年同月份的发病率。以2014年1—7月湖南省手足口病的实际月发病率对模型进行验证。表 3显示，预测值与实际值的相对误差绝对值范围为6.90%~46.31%，平均相对误差为20.37%。

表 3 (Table 3)

表 3 湖南省2014年1—7月手足口病月发病率预测值与实际值比较 日期 实际值(1/10万) 预测值(1/10万) 预测值95%CI(1/10万) 相对误差(%) 2014年1月 8.89 7.44 2.31~15.48 -16.33 2014年2月 5.02 5.48 0.68~14.91 9.31 2014年3月 19.68 10.88 2.64~24.74 -44.70 2014年4月 66.17 33.53 17.64~59.61 -46.31 2014年5月 59.57 65.73 39.68~98.31 10.34 2014年6月 55.61 50.75 27.96~80.27 -8.75 2014年7月 26.89 28.74 12.32~52.02 6.90

表 3 湖南省2014年1—7月手足口病月发病率预测值与实际值比较

图 6

Fig. 6

图 6 湖南省2008—2015年手足口病月发病率趋势预测图

表 4 (Table 4)

表 4 湖南省2014年8月—2015年7月手足口病月发病率预测值 日期 预测值(1/10万) 预测值95%CI(1/10万) 日期 预测值(1/10万) 预测值95%CI(1/10万) 2014年8月 18.52 5.97~38.00 2015年2月 7.76 0.70~22.41 2014年9月 24.11 9.27~45.92 2015年3月 13.49 2.92~31.79 2014年10月 26.92 11.04~49.86 2015年4月 31.70 13.37~57.83 2014年11月 27.88 11.62~51.14 2015年5月 46.74 23.59~77.72 2014年12月 24.53 9.49~46.58 2015年6月 38.60 17.90~67.17 2014年1月 11.01 1.95~27.45 2015年7月 26.32 9.90~50.62

表 4 湖南省2014年8月—2015年7月手足口病月发病率预测值

手足口病具有传染性强、传播途径复杂、传播速度快、隐性感染比例大、发病数多、无疫苗和特异性治疗方法等特点，是近年来备受关注的公共卫生问题。早期的预测和预警，可以了解手足口病疫情发病趋势，为防控提供重要的参考依据。

ARIMA模型只考虑时间序列的长期趋势，采取简单差分使得非平稳序列接近平稳；SARIMA模型在ARMIA模型的基础上，增加了对季节性和周期性的分析，考虑更加全面，可以采取季节性差分和(或)简单差分使得序列趋于平稳^{[11, 12]}。手足口病发病率时间序列往往具有长期趋势、周期性和季节性，适合用SARIMA模型拟合。SARIMA模型属于时间序列模型，基于历史数据来建模，不同省市各自建立的模型只适用于本省市，并不适用于其他地区。单次建立的SARIMA模型不能作为一成不变的预测工具^{[13, 14]}，只能用于短期预测。使用SARIMA模型长期地进行预测，需不断更新数据，重新拟合新的模型。在建立SARIMA模型的过程中，数据平稳化后各个备选模型的平稳R²等参数相差不大，提示数据的平稳化处理是关键步骤。不过，单从时间序列图、ACF以及PACF图来选择差分、季节性差分等处理方法，存在一定的主观性，故为了得到较好的模型，需要多次重复数据处理、模型参数评估等步骤^{[15, 16, 17]}。

本文采用SARIMA(1,0,0)(1,1,0)₁₂模型对湖南省手足口病月发病率趋势进行预测，预测值与实际值的平均相对误差为20.37%，低于黄晓霞等^[16]拟合的全国手足口病月发病率SARIMA(1,0,0)(0,1,0)₁₂模型的平均相对误差(28.62%)，拟合较好。通过月发病率趋势图可以直观、清楚得看出模型的拟合情况以及未来几个月湖南省手足口病疫情的发病情况、走向趋势。通过模型得到的预测值的95%CI可以用于预警，实际值若在95%CI以内，则疫情正常；若开始超过95%CI的上限(预警线)，则可以发出及时的预警，提醒疾病预防控制部门加强防控^[8]。本次研究的是单次建立的SARIMA模型，有较好的短期预测效果。为了得到较好的湖南省手足口病长期预测模型，有待后续研究中不断更新数据，重新拟合新的模型。