手足口病流行時間序列模型及其與氣象因素聯合預測研究*

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

【提 要】 目的 探討手足口病流行的時間序列特征與預測方法，為風險評估和政策措施制定提供科學依據。方法 收集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]。手足口病流行是一個復雜的自然和社會過程，氣象因素是手足口病感染的重要危險因素，非決定性因果關聯因素，但是可以借助氣象因素，更好地對手足口病流行做出預測。