手足口病流行时间序列模型及其与气象因素联合预测研究*

景钦隆 吴琦琳 鲁 影 陈 纯 张周斌

【提 要】 目的 探讨手足口病流行的时间序列特征与预测方法，为风险评估和政策措施制定提供科学依据。方法 收集2010至2017年广州市手足口病每月发病数和气象资料(平均气温、总和降雨量、相对湿度)，划分训练数据和验证数据，基于自回归求和移动平均模型(seasonal autoregressive integrated moving average，SARIMA)建立多元时间序列回归预测模型(SARIMA with external regressors，SARIMAX)。结果 广州市年平均发病人数为61795例，月平均发病人数为5150例。发病数时间序列具有明显的季节性特征，最终建立模型为SARIMA(0，1，0)(1，1，1)12纳入相对湿度滞后1期变量模型，该模型较基础SARIMA(0，1，0)(1，1，1)12模型AIC值降低16.52%，拟合过程均方根误差(root mean square error，RMSE)降低11.13%，预测过程RMSE降低40.68%。结论 SARIMAX模型可提高手足口病流行预测的精确度，相对湿度是广州地区手足口病流行的重要预测因素。

手足口病(hand-foot-mouth disease，HFMD)是5岁以下儿童多发传染病，大多数患者症状轻，以发热和手、足、口腔等部位的皮疹或疱疹为主要临床表现，重症患者病情凶险，病死风险高[1]。其病原肠道病毒血清型多，其中以肠道病毒71型和柯萨奇A16型为多见[2]。传播途径复杂，传染力强，季节性强，极易在幼托儿童聚集群体暴发流行。手足口病预测是日常防控工作的重要内容，与发病相关的危险因素较多，其中以气象因素最为关注[3]。

手足口病监测数据属于时间序列资料，因存在时间自相关性不满足独立性条件，传统线性回归等预测模型的适用性受限。季节自回归求和移动平均模型(seasonal autoregressive integrated moving average，SARIMA)因同时考虑序列自相关性、趋势和季节性得到广泛应用。近年来该模型进一步发展，基于该模型的多元时间序列回归模型(SARIMA with external regressors，SARIMAX)较单纯SARIMA模型显示了更好的预测效果，在传染病监测领域逐渐受到重视。本研究采用该模型方法对广州市手足口病发病人数时间序列结合气象因素联合预测，建立预测模型和识别重要气象因素，为手足口病预防控制提供科学依据。

资料与方法

1.资料来源及研究现场

选取2010年1月至2017年12月为研究时间段。每月手足口病发病人数来源于国家《传染病报告管理信息系统》。同期气象数据来源于中国气象数据共享服务网(http：//data.cma.cn/)，平均气温、总和降雨量和相对湿度按月度汇总。以2010年1月至2015年12月期间数据为模型拟合训练集，2016年1月至2017年12月数据为测试集。研究地点为广东省年发病人数最多的广州市[4]。

2.基础SARIMA乘积季节模型构建

SARIMA(p，d，q)(P，D，Q)S是一种广泛采用的时间序列分析和预测方法，其在原ARIMA模型基础上增加了对季节性和周期性的分析[5]。其中p为非季节自回归阶数，d为非季节差分阶数，q为非季节移动平均阶数，P为季节性自回归阶数，D为季节性差分阶数，Q为季节性移动平均阶数，S为季节性周期。

模型结构为[6]：

该模型综合考虑趋势、季节性和随机因素，对具有自相关性时间序列具有良好的预测效果[7]。建模过程如下：(1)序列平稳化：观察手足口病历年每月发病数时间序列图，如果具有明显的趋势和季节性，则通过序列差分和季节性差分方法，使序列成为零均值、无明显趋势且波动有界的平稳序列，纯随机性检验用LB(Ljung-Box)统计量。(2)模型识别与拟合：对差分后平稳序列，分别采用自相关函数(autocorrelation function，ACF)和偏自相关函数(partial autocorrelation function，PACF)初步确定p、d、q、P、D和Q的取值。S长度由实际专业背景分析得到。如果ACF结果拖尾，但PACF结果截尾，适合自回归(autoregressive model，AR)模型；如PACF拖尾，ACF截尾，则适合移动平均(moving average，MA)模型；如ACF和PACF均拖尾，则适合ARIMA模型。根据文献经验，P、Q取值不宜超过2阶，本研究取0、1和2逐个模型拟合。如果各参数取值困难，可列出可能的参数范围，多种参数组合模型逐一拟合，最后通过AIC(Akaike Information Criterion)值选择最佳模型，AIC越低，模型越佳[8]。本研究中模型参数采用最大似然估计(maximum likelihood estimation，MLE)法，模型系数检验采用t检验，显著性水平设定为α=0.05。

3.SARIMAX模型及预测分析

手足口病每月发病数与平均气温、总和降雨量和相对湿度相关关系存在滞后效应。本研究采用交叉相关函数计算与各气象因素不同滞后期相关系数，选择滞后相关系数最大且具有统计学意义(P<0.05)的滞后期变量纳入SARIMAX模型分析。本研究滞后期最大期数设定为4。

在选择最佳SARIMA基础模型后，将气象因素滞后变量(平均气温、总和降雨量、相对湿度)作为外部变量纳入该模型，即SARIMAX模型。模型结构如下[6]：

模型中X代表外部变量，可以纳入单因素变量，也可以纳入多因素变量，其他参数与基础模型结构参数意义一致。本次研究中，将通过自相关函数选择的滞后变量分别组合开展单因素和多因素SARIMAX分析。模型参数亦采用最大似然估计法估计，系数检验采用t检验，显著性水平为α=0.05。模型验证通过一步预测法比较模型预测的2016年1月至2017年12月的预测值与实际发病人数，计算均方根误差(root mean square error，RMSE)值。以模型系数均具有统计学意义，残差为纯随机序列，且预测误差RMSE低者为最佳模型[3]。

4.统计分析

本研究采用R语言软件(version3.4.4，the R foundation for Statistical Computing，Vienna，Austria)进行数据处理和模型分析，程序包应用包括tseries、forecast和TSA。

结 果

1.发病人数及气象因素概述

2010年至2017年，广州市累计报告手足口病494359例，年平均发病人数为61795例，月平均发病人数为(5150±4643)例。其中发病人数最多为2014年，共报告81152例(16.42%)，其次为2017年76559例(15.49%)、2013年72018例(14.57%)、2015年65218例(13.19%)、2016年60889例(12.32%)、2012年55284例(11.18%)、2011年46999例(9.51%)和2010年36240例(7.33%)。见图1。

平均气温、总和降雨量和相对湿度等气象因素亦具有明显的季节性，总体上数值表现为冬春季低，夏秋季高，季节性特征与发病特征类似。见图1。历年每月平均气温(21.97±5.62)℃，每月平均总和降雨量为(188.99±175.17)mm，每月平均相对湿度(77.76±6.82)%。

2.基础SARIMA乘积季节模型构建

图1a显示手足口病发病数时间序列在2011至2014年具有明显的趋势效应，且每年季节性效应明显。因此，在log转换后进行1阶12步差分得到近似平稳序列，LB统计量1至12阶P值均小于0.05，见图2。

根据差分变换的次数，初步确定模型为SARIMA(p，1，q)(P，1，Q)12，周期为12个月。ACF在0阶处截尾，同时在滞后12阶、24阶等周期性阶数有明显的波动，表明差分后序列仍然具有季节性效应。PACF在0阶处截尾，同时在滞后12阶处有明显的波动。因此考虑p的取值为0或1，q的取值为0或1。P、Q的取值判断，参考既往文献经验不超过2，因此取值0、1、2分别逐个不同组合模型测试。合计36个不同组合模型中，10个模型系数有统计学意义，其中SARIMA(0，1，0)(1，1，1)12模型的AIC值(78.78)最小，为拟合最佳模型，决定系数R2为90.25%。见表1。

图1 手足口病历年每月发病数及气象因素时间序列图a)历年每月发病数序列图，b)历年每月总和降雨量序列图，c)历年每月平均气温序列图，d)历年每月平均相对湿度序列图

图2 手足口病每月发病数log转换后1阶12步差分序列及其ACF和PACF图a)每月发病数log转换后1阶12步差分序列图，b)发病数1阶12步差分后自相关图，c)发病数1阶12步差分后偏自相关图

模型010111010012010210010110111212011011010011011010110010010010AIC78.779780.467181.063482.313882.826383.266284.1218117.6920117.9636119.1907

*：模型一行中6位数字分别代表p、d、q、P、D和Q值。

3.SARIMAX模型与预测分析

以每月发病数(log转换)为目标变量，交叉相关分析显示平均气温滞后0期的相关系数最大(r=0.7930，P<0.05)，总和降雨量滞后1期的相关系数最大(r=0.6058，P<0.05)，平均相对湿度滞后1期的相关系数最大(r=0.7468，P<0.05)。见表2。

表2 手足口病每月发病数(log转换)与气象因素交叉相关分析表(α=0.05)

*：“+”为正相关且具有统计学意义，“.”没有统计学意义，“-”负相关且具有统计学意义

将交叉相关分析结果中，选择与平均气温滞后0期、总和降雨量滞后1期和平均相对湿度滞后1期交叉相关系数最大且具有统计学意义的滞后期变量纳入SARIMAX模型分析。单因素和多因素分析结果显示，SARIMA(0，1，0)(1，1，1)12模型纳入平均相对湿度滞后1期的单因素模型拟合AIC值最小(65.7638)，决定系数R2为90.65%，且模型系数均具有统计学意义(P<0.05)，其它单因素和多因素模型拟合模型系数在α=0.05显著性水平上不全有统计学意义。见表3。SARIMA(0，1，0)(1，1，1)12纳入平均相对湿度滞后1期变量模型较基础SARIMA(0，1，0)(1，1，1)12模型AIC值降低16.52%，拟合过程RMSE降低11.13%，预测过程RMSE降低40.68%。模型残差均为随机序列，LB统计量1至12阶P值均大于0.05。

表3 手足口病每月发病数(log转换)与气象因素ARIMAX模型分析结果

*:sar1：季节性AR(1)，sma1：季节性MA(2)，“**”：P< 0.01，“*”：P<0.05，T(lag0)：平均气温滞后0期，P(lag1)：总和降雨量滞后1期，RH(lag1)：平均相对湿度滞后1期。②model0为基础模型，SARIMA(0，1，0)(1，1，1)12模型；model1至model7分别为基础模型纳入：model1，纳入T(lag0)；model2，纳入P(lag1)；model3，纳入RH(lag1)；model4，纳入T(lag0)、P(lag1)；model5，纳入T(lag0)、RH(lag1)；model6，纳入P(lag1)、RH(lag1)；model7，纳入T(lag0)、P(lag1)、RH(lag1)。

讨 论

本研究建立的基于SARIMA模型的多元时间序列回归模型SARIMA(0，1，0)(1，1，1)12纳入平均相对湿度滞后1期模型较单纯SARIMA(0，1，0)(1，1，1)12模型在拟合优度和预测效果方面均具有更好的精确度，特别是在短期预测方面[9]，这在一定程度上克服了后者未考虑外部变量因素造成的准确性和可扩展性问题[8]。该模型可为手足口病风险评估和政策措施制定提供简单有效的手段，有助于提前应对和降低流行强度。

我国南北方在手足口病流行季节性方面有所差异，北方地区峰值在6月份，南方地区则表现为大小两个高峰[2，10]，广州市手足口病流行亦表现为大小2个高峰[11]。目前认为，手足口病流行季节性周期原因与气象因素关系密切。既往研究显示气温、降雨量、相对湿度、气压、风速、日照等因素与该病均有不同程度相关性，且不同研究者结果存在差异[12]。本研究，平均气温和降雨量纳入模型后，增大了拟合优度值AIC，而且模型参数系数亦不具有统计学意义，最终仅将相对湿度滞后1期纳入模型，达到了模型简单，预测效果也得以提升的目的。

SARIMAX同时是识别传染病危险因素的重要方法。郑州手足口病流行的研究中，气温滞后2周是手足口病发病的关键危险因素，其它气象因素如相对湿度、日照时间、降雨量和风速纳入模型则无统计学意义[12-13]。本研究发现相对湿度是手足口病发病的关键影响因素，而温度和降雨量则无统计学意义。相对湿度在同为肠道传染病的脊髓灰质炎流行中是关键气象因素[14]。广州市3月至6月相对湿度处于全年较高水平，该时段滞后1月手足口病发病数亦居年内最高位。原因可能与较高的相对湿度，有利于肠道病毒延长生存时间和有利于增殖过程[15]有关。其次，该时段为幼托儿童在校期，人口在幼托场所相对密集，传播途径容易实现，增加了病毒感染的机会。在同样为在校期的9月至11月，因相对湿度处于较低水平，此时间段的流行峰值明显低于4月至7月。

本次研究亦存在如下局限性：本次研究外部变量因素仅纳入气象因素，实际上影响手足口病的的因素众多，比如人口经济学因素、人口密度与流动、防控政策与措施落实、病原血清型别转换等问题；其次，国家传染病报告管理系统，因一线医生的诊断报告意识、轻症患者不就医和隐性感染者，存在不能识别整个疾病谱问题；该模型在短期预测方面效果较好，随着新的数据不断产生，模型需要一定的调整才能适应预测需要[16]。手足口病流行是一个复杂的自然和社会过程，气象因素是手足口病感染的重要危险因素，非决定性因果关联因素，但是可以借助气象因素，更好地对手足口病流行做出预测。