手足口病流行時(shí)間序列模型及其與氣象因素聯(lián)合預(yù)測(cè)研究*

景欽隆 吳琦琳 魯 影 陳 純 張周斌

【提 要】 目的 探討手足口病流行的時(shí)間序列特征與預(yù)測(cè)方法，為風(fēng)險(xiǎn)評(píng)估和政策措施制定提供科學(xué)依據(jù)。方法 收集2010至2017年廣州市手足口病每月發(fā)病數(shù)和氣象資料(平均氣溫、總和降雨量、相對(duì)濕度)，劃分訓(xùn)練數(shù)據(jù)和驗(yàn)證數(shù)據(jù)，基于自回歸求和移動(dòng)平均模型(seasonal autoregressive integrated moving average，SARIMA)建立多元時(shí)間序列回歸預(yù)測(cè)模型(SARIMA with external regressors，SARIMAX)。結(jié)果 廣州市年平均發(fā)病人數(shù)為61795例，月平均發(fā)病人數(shù)為5150例。發(fā)病數(shù)時(shí)間序列具有明顯的季節(jié)性特征，最終建立模型為SARIMA(0，1，0)(1，1，1)12納入相對(duì)濕度滯后1期變量模型，該模型較基礎(chǔ)SARIMA(0，1，0)(1，1，1)12模型AIC值降低16.52%，擬合過(guò)程均方根誤差(root mean square error，RMSE)降低11.13%，預(yù)測(cè)過(guò)程RMSE降低40.68%。結(jié)論 SARIMAX模型可提高手足口病流行預(yù)測(cè)的精確度，相對(duì)濕度是廣州地區(qū)手足口病流行的重要預(yù)測(cè)因素。

手足口病(hand-foot-mouth disease，HFMD)是5歲以下兒童多發(fā)傳染病，大多數(shù)患者癥狀輕，以發(fā)熱和手、足、口腔等部位的皮疹或皰疹為主要臨床表現(xiàn)，重癥患者病情兇險(xiǎn)，病死風(fēng)險(xiǎn)高[1]。其病原腸道病毒血清型多，其中以腸道病毒71型和柯薩奇A16型為多見(jiàn)[2]。傳播途徑復(fù)雜，傳染力強(qiáng)，季節(jié)性強(qiáng)，極易在幼托兒童聚集群體暴發(fā)流行。手足口病預(yù)測(cè)是日常防控工作的重要內(nèi)容，與發(fā)病相關(guān)的危險(xiǎn)因素較多，其中以氣象因素最為關(guān)注[3]。

手足口病監(jiān)測(cè)數(shù)據(jù)屬于時(shí)間序列資料，因存在時(shí)間自相關(guān)性不滿(mǎn)足獨(dú)立性條件，傳統(tǒng)線(xiàn)性回歸等預(yù)測(cè)模型的適用性受限。季節(jié)自回歸求和移動(dòng)平均模型(seasonal autoregressive integrated moving average，SARIMA)因同時(shí)考慮序列自相關(guān)性、趨勢(shì)和季節(jié)性得到廣泛應(yīng)用。近年來(lái)該模型進(jìn)一步發(fā)展，基于該模型的多元時(shí)間序列回歸模型(SARIMA with external regressors，SARIMAX)較單純SARIMA模型顯示了更好的預(yù)測(cè)效果，在傳染病監(jiān)測(cè)領(lǐng)域逐漸受到重視。本研究采用該模型方法對(duì)廣州市手足口病發(fā)病人數(shù)時(shí)間序列結(jié)合氣象因素聯(lián)合預(yù)測(cè)，建立預(yù)測(cè)模型和識(shí)別重要?dú)庀笠蛩?，為手足口病預(yù)防控制提供科學(xué)依據(jù)。

資料與方法

1.資料來(lái)源及研究現(xiàn)場(chǎng)

選取2010年1月至2017年12月為研究時(shí)間段。每月手足口病發(fā)病人數(shù)來(lái)源于國(guó)家《傳染病報(bào)告管理信息系統(tǒng)》。同期氣象數(shù)據(jù)來(lái)源于中國(guó)氣象數(shù)據(jù)共享服務(wù)網(wǎng)(http：//data.cma.cn/)，平均氣溫、總和降雨量和相對(duì)濕度按月度匯總。以2010年1月至2015年12月期間數(shù)據(jù)為模型擬合訓(xùn)練集，2016年1月至2017年12月數(shù)據(jù)為測(cè)試集。研究地點(diǎn)為廣東省年發(fā)病人數(shù)最多的廣州市[4]。

2.基礎(chǔ)SARIMA乘積季節(jié)模型構(gòu)建

SARIMA(p，d，q)(P，D，Q)S是一種廣泛采用的時(shí)間序列分析和預(yù)測(cè)方法，其在原ARIMA模型基礎(chǔ)上增加了對(duì)季節(jié)性和周期性的分析[5]。其中p為非季節(jié)自回歸階數(shù)，d為非季節(jié)差分階數(shù)，q為非季節(jié)移動(dòng)平均階數(shù)，P為季節(jié)性自回歸階數(shù)，D為季節(jié)性差分階數(shù)，Q為季節(jié)性移動(dòng)平均階數(shù)，S為季節(jié)性周期。

模型結(jié)構(gòu)為[6]：

該模型綜合考慮趨勢(shì)、季節(jié)性和隨機(jī)因素，對(duì)具有自相關(guān)性時(shí)間序列具有良好的預(yù)測(cè)效果[7]。建模過(guò)程如下：(1)序列平穩(wěn)化：觀(guān)察手足口病歷年每月發(fā)病數(shù)時(shí)間序列圖，如果具有明顯的趨勢(shì)和季節(jié)性，則通過(guò)序列差分和季節(jié)性差分方法，使序列成為零均值、無(wú)明顯趨勢(shì)且波動(dòng)有界的平穩(wěn)序列，純隨機(jī)性檢驗(yàn)用LB(Ljung-Box)統(tǒng)計(jì)量。(2)模型識(shí)別與擬合：對(duì)差分后平穩(wěn)序列，分別采用自相關(guān)函數(shù)(autocorrelation function，ACF)和偏自相關(guān)函數(shù)(partial autocorrelation function，PACF)初步確定p、d、q、P、D和Q的取值。S長(zhǎng)度由實(shí)際專(zhuān)業(yè)背景分析得到。如果ACF結(jié)果拖尾，但PACF結(jié)果截尾，適合自回歸(autoregressive model，AR)模型；如PACF拖尾，ACF截尾，則適合移動(dòng)平均(moving average，MA)模型；如ACF和PACF均拖尾，則適合ARIMA模型。根據(jù)文獻(xiàn)經(jīng)驗(yàn)，P、Q取值不宜超過(guò)2階，本研究取0、1和2逐個(gè)模型擬合。如果各參數(shù)取值困難，可列出可能的參數(shù)范圍，多種參數(shù)組合模型逐一擬合，最后通過(guò)AIC(Akaike Information Criterion)值選擇最佳模型，AIC越低，模型越佳[8]。本研究中模型參數(shù)采用最大似然估計(jì)(maximum likelihood estimation，MLE)法，模型系數(shù)檢驗(yàn)采用t檢驗(yàn)，顯著性水平設(shè)定為α=0.05。

3.SARIMAX模型及預(yù)測(cè)分析

手足口病每月發(fā)病數(shù)與平均氣溫、總和降雨量和相對(duì)濕度相關(guān)關(guān)系存在滯后效應(yīng)。本研究采用交叉相關(guān)函數(shù)計(jì)算與各氣象因素不同滯后期相關(guān)系數(shù)，選擇滯后相關(guān)系數(shù)最大且具有統(tǒng)計(jì)學(xué)意義(P<0.05)的滯后期變量納入SARIMAX模型分析。本研究滯后期最大期數(shù)設(shè)定為4。

在選擇最佳SARIMA基礎(chǔ)模型后，將氣象因素滯后變量(平均氣溫、總和降雨量、相對(duì)濕度)作為外部變量納入該模型，即SARIMAX模型。模型結(jié)構(gòu)如下[6]：

模型中X代表外部變量，可以納入單因素變量，也可以納入多因素變量，其他參數(shù)與基礎(chǔ)模型結(jié)構(gòu)參數(shù)意義一致。本次研究中，將通過(guò)自相關(guān)函數(shù)選擇的滯后變量分別組合開(kāi)展單因素和多因素SARIMAX分析。模型參數(shù)亦采用最大似然估計(jì)法估計(jì)，系數(shù)檢驗(yàn)采用t檢驗(yàn)，顯著性水平為α=0.05。模型驗(yàn)證通過(guò)一步預(yù)測(cè)法比較模型預(yù)測(cè)的2016年1月至2017年12月的預(yù)測(cè)值與實(shí)際發(fā)病人數(shù)，計(jì)算均方根誤差(root mean square error，RMSE)值。以模型系數(shù)均具有統(tǒng)計(jì)學(xué)意義，殘差為純隨機(jī)序列，且預(yù)測(cè)誤差RMSE低者為最佳模型[3]。

4.統(tǒng)計(jì)分析

本研究采用R語(yǔ)言軟件(version3.4.4，the R foundation for Statistical Computing，Vienna，Austria)進(jìn)行數(shù)據(jù)處理和模型分析，程序包應(yīng)用包括tseries、forecast和TSA。

結(jié) 果

1.發(fā)病人數(shù)及氣象因素概述

2010年至2017年，廣州市累計(jì)報(bào)告手足口病494359例，年平均發(fā)病人數(shù)為61795例，月平均發(fā)病人數(shù)為(5150±4643)例。其中發(fā)病人數(shù)最多為2014年，共報(bào)告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%)。見(jiàn)圖1。

平均氣溫、總和降雨量和相對(duì)濕度等氣象因素亦具有明顯的季節(jié)性，總體上數(shù)值表現(xiàn)為冬春季低，夏秋季高，季節(jié)性特征與發(fā)病特征類(lèi)似。見(jiàn)圖1。歷年每月平均氣溫(21.97±5.62)℃，每月平均總和降雨量為(188.99±175.17)mm，每月平均相對(duì)濕度(77.76±6.82)%。

2.基礎(chǔ)SARIMA乘積季節(jié)模型構(gòu)建

圖1a顯示手足口病發(fā)病數(shù)時(shí)間序列在2011至2014年具有明顯的趨勢(shì)效應(yīng)，且每年季節(jié)性效應(yīng)明顯。因此，在log轉(zhuǎn)換后進(jìn)行1階12步差分得到近似平穩(wěn)序列，LB統(tǒng)計(jì)量1至12階P值均小于0.05，見(jiàn)圖2。

根據(jù)差分變換的次數(shù)，初步確定模型為SARIMA(p，1，q)(P，1，Q)12，周期為12個(gè)月。ACF在0階處截尾，同時(shí)在滯后12階、24階等周期性階數(shù)有明顯的波動(dòng)，表明差分后序列仍然具有季節(jié)性效應(yīng)。PACF在0階處截尾，同時(shí)在滯后12階處有明顯的波動(dòng)。因此考慮p的取值為0或1，q的取值為0或1。P、Q的取值判斷，參考既往文獻(xiàn)經(jīng)驗(yàn)不超過(guò)2，因此取值0、1、2分別逐個(gè)不同組合模型測(cè)試。合計(jì)36個(gè)不同組合模型中，10個(gè)模型系數(shù)有統(tǒng)計(jì)學(xué)意義，其中SARIMA(0，1，0)(1，1，1)12模型的AIC值(78.78)最小，為擬合最佳模型，決定系數(shù)R2為90.25%。見(jiàn)表1。

圖1 手足口病歷年每月發(fā)病數(shù)及氣象因素時(shí)間序列圖a)歷年每月發(fā)病數(shù)序列圖，b)歷年每月總和降雨量序列圖，c)歷年每月平均氣溫序列圖，d)歷年每月平均相對(duì)濕度序列圖

圖2 手足口病每月發(fā)病數(shù)log轉(zhuǎn)換后1階12步差分序列及其ACF和PACF圖a)每月發(fā)病數(shù)log轉(zhuǎn)換后1階12步差分序列圖，b)發(fā)病數(shù)1階12步差分后自相關(guān)圖，c)發(fā)病數(shù)1階12步差分后偏自相關(guān)圖

模型010111010012010210010110111212011011010011011010110010010010AIC78.779780.467181.063482.313882.826383.266284.1218117.6920117.9636119.1907

*：模型一行中6位數(shù)字分別代表p、d、q、P、D和Q值。

3.SARIMAX模型與預(yù)測(cè)分析

以每月發(fā)病數(shù)(log轉(zhuǎn)換)為目標(biāo)變量，交叉相關(guān)分析顯示平均氣溫滯后0期的相關(guān)系數(shù)最大(r=0.7930，P<0.05)，總和降雨量滯后1期的相關(guān)系數(shù)最大(r=0.6058，P<0.05)，平均相對(duì)濕度滯后1期的相關(guān)系數(shù)最大(r=0.7468，P<0.05)。見(jiàn)表2。

表2 手足口病每月發(fā)病數(shù)(log轉(zhuǎn)換)與氣象因素交叉相關(guān)分析表(α=0.05)

*：“+”為正相關(guān)且具有統(tǒng)計(jì)學(xué)意義，“.”沒(méi)有統(tǒng)計(jì)學(xué)意義，“-”負(fù)相關(guān)且具有統(tǒng)計(jì)學(xué)意義

將交叉相關(guān)分析結(jié)果中，選擇與平均氣溫滯后0期、總和降雨量滯后1期和平均相對(duì)濕度滯后1期交叉相關(guān)系數(shù)最大且具有統(tǒng)計(jì)學(xué)意義的滯后期變量納入SARIMAX模型分析。單因素和多因素分析結(jié)果顯示，SARIMA(0，1，0)(1，1，1)12模型納入平均相對(duì)濕度滯后1期的單因素模型擬合AIC值最小(65.7638)，決定系數(shù)R2為90.65%，且模型系數(shù)均具有統(tǒng)計(jì)學(xué)意義(P<0.05)，其它單因素和多因素模型擬合模型系數(shù)在α=0.05顯著性水平上不全有統(tǒng)計(jì)學(xué)意義。見(jiàn)表3。SARIMA(0，1，0)(1，1，1)12納入平均相對(duì)濕度滯后1期變量模型較基礎(chǔ)SARIMA(0，1，0)(1，1，1)12模型AIC值降低16.52%，擬合過(guò)程RMSE降低11.13%，預(yù)測(cè)過(guò)程RMSE降低40.68%。模型殘差均為隨機(jī)序列，LB統(tǒng)計(jì)量1至12階P值均大于0.05。

表3 手足口病每月發(fā)病數(shù)(log轉(zhuǎn)換)與氣象因素ARIMAX模型分析結(jié)果

*:sar1：季節(jié)性AR(1)，sma1：季節(jié)性MA(2)，“**”：P< 0.01，“*”：P<0.05，T(lag0)：平均氣溫滯后0期，P(lag1)：總和降雨量滯后1期，RH(lag1)：平均相對(duì)濕度滯后1期。②model0為基礎(chǔ)模型，SARIMA(0，1，0)(1，1，1)12模型；model1至model7分別為基礎(chǔ)模型納入：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模型的多元時(shí)間序列回歸模型SARIMA(0，1，0)(1，1，1)12納入平均相對(duì)濕度滯后1期模型較單純SARIMA(0，1，0)(1，1，1)12模型在擬合優(yōu)度和預(yù)測(cè)效果方面均具有更好的精確度，特別是在短期預(yù)測(cè)方面[9]，這在一定程度上克服了后者未考慮外部變量因素造成的準(zhǔn)確性和可擴(kuò)展性問(wèn)題[8]。該模型可為手足口病風(fēng)險(xiǎn)評(píng)估和政策措施制定提供簡(jiǎn)單有效的手段，有助于提前應(yīng)對(duì)和降低流行強(qiáng)度。

我國(guó)南北方在手足口病流行季節(jié)性方面有所差異，北方地區(qū)峰值在6月份，南方地區(qū)則表現(xiàn)為大小兩個(gè)高峰[2，10]，廣州市手足口病流行亦表現(xiàn)為大小2個(gè)高峰[11]。目前認(rèn)為，手足口病流行季節(jié)性周期原因與氣象因素關(guān)系密切。既往研究顯示氣溫、降雨量、相對(duì)濕度、氣壓、風(fēng)速、日照等因素與該病均有不同程度相關(guān)性，且不同研究者結(jié)果存在差異[12]。本研究，平均氣溫和降雨量納入模型后，增大了擬合優(yōu)度值A(chǔ)IC，而且模型參數(shù)系數(shù)亦不具有統(tǒng)計(jì)學(xué)意義，最終僅將相對(duì)濕度滯后1期納入模型，達(dá)到了模型簡(jiǎn)單，預(yù)測(cè)效果也得以提升的目的。

SARIMAX同時(shí)是識(shí)別傳染病危險(xiǎn)因素的重要方法。鄭州手足口病流行的研究中，氣溫滯后2周是手足口病發(fā)病的關(guān)鍵危險(xiǎn)因素，其它氣象因素如相對(duì)濕度、日照時(shí)間、降雨量和風(fēng)速納入模型則無(wú)統(tǒng)計(jì)學(xué)意義[12-13]。本研究發(fā)現(xiàn)相對(duì)濕度是手足口病發(fā)病的關(guān)鍵影響因素，而溫度和降雨量則無(wú)統(tǒng)計(jì)學(xué)意義。相對(duì)濕度在同為腸道傳染病的脊髓灰質(zhì)炎流行中是關(guān)鍵氣象因素[14]。廣州市3月至6月相對(duì)濕度處于全年較高水平，該時(shí)段滯后1月手足口病發(fā)病數(shù)亦居年內(nèi)最高位。原因可能與較高的相對(duì)濕度，有利于腸道病毒延長(zhǎng)生存時(shí)間和有利于增殖過(guò)程[15]有關(guān)。其次，該時(shí)段為幼托兒童在校期，人口在幼托場(chǎng)所相對(duì)密集，傳播途徑容易實(shí)現(xiàn)，增加了病毒感染的機(jī)會(huì)。在同樣為在校期的9月至11月，因相對(duì)濕度處于較低水平，此時(shí)間段的流行峰值明顯低于4月至7月。

本次研究亦存在如下局限性：本次研究外部變量因素僅納入氣象因素，實(shí)際上影響手足口病的的因素眾多，比如人口經(jīng)濟(jì)學(xué)因素、人口密度與流動(dòng)、防控政策與措施落實(shí)、病原血清型別轉(zhuǎn)換等問(wèn)題；其次，國(guó)家傳染病報(bào)告管理系統(tǒng)，因一線(xiàn)醫(yī)生的診斷報(bào)告意識(shí)、輕癥患者不就醫(yī)和隱性感染者，存在不能識(shí)別整個(gè)疾病譜問(wèn)題；該模型在短期預(yù)測(cè)方面效果較好，隨著新的數(shù)據(jù)不斷產(chǎn)生，模型需要一定的調(diào)整才能適應(yīng)預(yù)測(cè)需要[16]。手足口病流行是一個(gè)復(fù)雜的自然和社會(huì)過(guò)程，氣象因素是手足口病感染的重要危險(xiǎn)因素，非決定性因果關(guān)聯(lián)因素，但是可以借助氣象因素，更好地對(duì)手足口病流行做出預(yù)測(cè)。