應(yīng)用SARIMA和ETS模型對湖南省腎綜合征出血熱發(fā)病趨勢的預(yù)測*

哈爾濱醫(yī)科大學(xué)衛(wèi)生管理學(xué)院社會醫(yī)學(xué)與衛(wèi)生事業(yè)管理教研室(150081) 邵瑛琦 劉 歡 李晨希 孟祥偉 李 樂 王 星 吳群紅

【提 要】 目的 探討比較SARIMA模型和ETS模型在湖南省腎綜合征出血熱的發(fā)病預(yù)測的應(yīng)用，為腎綜合征出血熱的防控提供依據(jù)。方法 利用2005-2014年腎綜合征出血熱月度發(fā)病數(shù)據(jù)建立SARIMA模型和ETS模型，并通過模型預(yù)測2015年1～12月的腎綜合征出血熱發(fā)病數(shù)，用均方根誤差(root mean square error，RMSE)和絕對百分比誤差(absolutc percentage error，MAPE)作為評價指標。結(jié)果 SARIMA(1，0，0)(3，0，0)12模型是腎綜合征出血熱發(fā)病趨勢的最優(yōu)擬合預(yù)測模型，SARIMA模型的MAPE為13.18%，低于ETS模型31.14%，SARIMA模型的RMSE為16.27%也低于ETS模型的25.88%。結(jié)論 SARIMA(1，0，0)(3，0，0)12模型模擬擬合效果較好，預(yù)測結(jié)果可為今后腎綜合征出血熱的預(yù)防和控制提供理論支持。

腎綜合征出血熱(hemorrhagic fever with renal syndrome，HFRS)是由漢坦病毒(hantaviruses，HV)引起的一種自然疫源性疾病，鼠類為自然宿主和主要傳染源，臨床表現(xiàn)以高熱、腰痛、眼眶痛及腎功能損害為特征[1]。我國是全球HFRS發(fā)病最多的國家，累計報告病例占全球病例總數(shù)的90%以上[2]。2014年監(jiān)測結(jié)果顯示，湖南省為腎綜合征出血熱高發(fā)省份之一，全年發(fā)病584人[3]。為了解湖南省HFRS流行特點及未來的流行趨勢，本研究通過季節(jié)自回歸移動平均(seasonal auto regressive integrated moving average，SARIMA)模型和指數(shù)平滑模型(exponential smoothing，ETS)對湖南省 HFRS流行趨勢進行預(yù)測，探討 SARIMA模型和ETS模型在 HFRS中的預(yù)測效果，為湖南省HFRS防控提供科學(xué)依據(jù)[4]。

資料與方法

1.資料來源

資料來源于公共衛(wèi)生科學(xué)數(shù)據(jù)中心(http：//www.phsciencedata.cn/)中湖南省HFRS 2005-2015年的月度數(shù)據(jù)，根據(jù)逐月發(fā)病數(shù)據(jù)建立模型，2015年1～12月HFRS逐月發(fā)病人數(shù)實際值驗證模型預(yù)測效果。

2.方法

使用Excel 2016建立HFRS數(shù)據(jù)庫，利用R 3.4.2軟件。利用R軟件對湖南省2005-2015年的HFRS月度發(fā)病數(shù)據(jù)進行處理和分析。選擇2005年1月-2015年12月的數(shù)據(jù)進行建模，通過指數(shù)平滑模型ETS和SARIMA模型對2015年發(fā)病人數(shù)進行預(yù)測。

(1)指數(shù)平滑模型(ETS)

基于Hyndman等人描述的方法及R軟件環(huán)境中的預(yù)測包建立ETS預(yù)測模型。指數(shù)平滑法是用序列以往實際觀測值的加權(quán)平均來預(yù)測未來數(shù)值，序列中近期的數(shù)據(jù)賦予較大的權(quán)重，遠期的數(shù)據(jù)賦予較小的權(quán)重[5]。指數(shù)平滑法通過誤差(error)、整體趨勢(trend)和季節(jié)性(seasonal)三個主要參數(shù)之間的相加、相乘或無運算進行模型擬合。在參數(shù)選擇過程中，自動選擇ETS模型來擬合具有可乘成分的指數(shù)模型，并評估替代模型以選擇性能最佳的模型模擬數(shù)據(jù)。在模型選擇過程中，基于赤遲信息準則(Akaike information criterion，AIC)的最小值、校正的Akaike信息準則(the corrected Akaike information criterion，AICc)或貝葉斯信息準則(Bayesian information criterion，BIC)選擇最佳模型。Ljung-Box Q檢驗用于診斷殘留錯誤序列是否為白噪聲序列。

(2)季節(jié)性自回歸綜合移動平均模型(SARIMA)

季節(jié)性差分自回歸滑動平均模型(seasonal autoregressive integrated moving average，SARIMA)建立在差分自回歸移動平均模型(autoregressive integrated moving average model，ARIMA)[6]基礎(chǔ)上，增加了對季節(jié)性和周期性的分析，表達為SARIMA(p，d，q)′(P，D，Q)s，其中p為自回歸階數(shù)，d為一般差分階數(shù)，q為移動平均階數(shù)，P為季節(jié)性自回歸階數(shù)，D為季節(jié)性差分階數(shù)，Q為季節(jié)性移動平均階數(shù)，s為季節(jié)周期[7]。SARIMA(p，q)模型的數(shù)學(xué)表達式為：Yt=μ+β1Yt-1+β2Yt-2+…+βpYt-p-θ1εt-1-θ2εt-2，-…-θqεq-2+εt。其中Yt是級數(shù)的給定值，β和θ是權(quán)重，μ是級數(shù)的平均值，而εt是預(yù)測誤差。其中k表示模型中將包括的最大時滯數(shù)，而Yt是經(jīng)過差分后獲得的新序列，φ是相應(yīng)的系數(shù)。

本研究使用R 3.4.2軟件進行SARIMA 模型的建模，其過程主要包括序列平穩(wěn)性檢驗、模型識別、模型的參數(shù)估計、模型的檢驗和優(yōu)化。①序列樣本平穩(wěn)化檢驗：通過差分使非平穩(wěn)時間序列平穩(wěn)。預(yù)測包中的ndiffs函數(shù)可用于幫助確定d的最佳值。利用單位根(augmented dickey-fuller，ADF)檢驗進行平穩(wěn)化處理。當ADF測試的結(jié)果為P<0.05，則證明該序列是穩(wěn)定的。②模型識別：根據(jù)自相關(guān)函數(shù)(autocorrelation function，ACF)圖和偏自相關(guān)(partial autocorrelation，PACF)圖估算SARIMA模型的參數(shù)。③參數(shù)估計和檢驗：利用預(yù)測包中的auto.arima函數(shù)根據(jù)AIC、AICc或BIC中的最小值來選擇最佳的SARIMA模型。最后還要對模型的參數(shù)進行顯著性檢驗和對殘差進行白噪聲檢驗，使用Ljung-Box Q檢驗來驗證估計的殘差是否滿足白噪聲序列的要求，若殘差為非白噪聲序列，則需要對模型的參數(shù)重新估計和擬合。

(3)效果評估指標

采用均方根誤差(root mean squared error，RMSE)和絕對百分比誤差(mean absolute percentage error，MAPE)2個指標評價ETS和SARIMA模型的擬合、預(yù)測效果，其中，2個指標均最小者，擬合(預(yù)測)效果最好。

結(jié) 果

1.2005年至2015年HFRS報告發(fā)病數(shù)和時空分布情況

湖南省2005-2015年HFRS月發(fā)病情況見表1?；诒?的數(shù)據(jù)運用R語言中的“decompose”函數(shù)分解時間序列的趨勢性、季節(jié)性和隨機誤差[8]并繪制時間序列分解圖(圖1)。時序圖顯示了2005年1月至2014年12月期間每月的HFRS發(fā)病趨勢。HFRS發(fā)病情況季節(jié)特征明顯，呈雙峰型。春季為HFRS的第一個發(fā)病高峰期，隨后有短暫的下降趨勢，秋季面臨一個發(fā)病小高峰。2005-2010年與2011-2014年發(fā)病人數(shù)均呈現(xiàn)出先升高后降低的趨勢，發(fā)病最高峰為2013年(圖1)。

表1 2005-2015年湖南省HFRS逐月發(fā)病人數(shù)

圖1 2005年1月至2014年12月湖南省HFRS發(fā)病時間序列分解圖

2.模型檢驗與處理

軟件中的nddifs代碼結(jié)果顯示HFRS時間序列不需要差分平穩(wěn)化處理(d=0)。ADF測試結(jié)果表明統(tǒng)計學(xué)上有意義(P=0.01)，這表明時間序列是穩(wěn)定的。對模型進行殘差檢驗以估算其他參數(shù)，繪制殘差序列的自相關(guān)(ACF)和偏自相關(guān)圖(PACF)(圖2)。在樣本數(shù)據(jù)中使用時間序列建模器構(gòu)建SARIMA模型，軟件自動選擇SARIMA(0，0，1)(3，0，0)12作為最佳擬合模型(AIC=1025.19，BIC=1041.92，AICc=1025.94)。Ljung-Box Q檢驗進一步表明，殘留誤差序列包含白噪聲(χ2=0.014229，df=1，P=0.905)。在運行ETS代碼時，軟件自動選擇ETS(M，N，M)模型(AIC=1199.165，BIC=1240.977，AICc=1203.780)作為性能最佳的ETS模型。ETS(M，N，M)模型的Ljung-Box Q測試結(jié)果表明，殘留誤差序列也包含白噪聲(χ2=0.45363，df=1，P=0.5006)。

圖2 差分HFRS時間序列的自相關(guān)和偏自相關(guān)圖

3.模型預(yù)測

運用模型SARIMA(0，0，1)(3，0，0)12和ETS(M，N，M)模型對2015年1月至2015年12月的HFRS疾病逐月發(fā)病人數(shù)進行回代預(yù)測，得到發(fā)病人數(shù)預(yù)測表(表2)和發(fā)病趨勢圖(圖3)，結(jié)果顯示SARIMA模型和EST模型擬合的2015年每月HFRS實際觀測發(fā)病人數(shù)均落入95%的置信區(qū)間之內(nèi)。

圖3 2015年湖南省HFRS發(fā)病趨勢預(yù)測

表2 湖南省2015年1～12月HFRS發(fā)病人數(shù)預(yù)測

4.評價擬合效果

SARIMA(1，0，0)(3，0，0)12和ETS(M，N，M)模型的性能度量結(jié)果如表3所示。從仿真的兩個方面比較兩個模型的性能結(jié)果表明，SARIMA(1，0，0)(3，0，0)12模型的均方根誤差(RMSE)和平均絕對百分比誤差(MAPE)均低于ETS(M，N，M)模型?；诤鲜?005年1月-2015年12月短期HFRS預(yù)測結(jié)果，表明SARIMA模型的預(yù)測值與真實值更為吻合，模擬效果較好。

表3 SARIMAM模型、ETS模型與真實值對比

討 論

傳染病具有一定的發(fā)病規(guī)律，如線性趨勢、季節(jié)性、周期性等，通過時間序列分析數(shù)據(jù)探索傳染病在時間上的變化規(guī)律，能夠預(yù)測其未來發(fā)展趨勢[9]。總體來看湖南省HFRS在2005-2015年期間發(fā)病表現(xiàn)為上升到下降，再上升至下降的波動狀態(tài)，HFRS的防控取得了一定的成效，疫情有所緩解，這和相關(guān)報道一致[10]。但因存在疫情反復(fù)的可能，仍是當前需要重點關(guān)注的公共衛(wèi)生問題[11]。

隨著預(yù)測理論及預(yù)測技術(shù)的發(fā)展與完善，越來越多的模型被應(yīng)用于探索傳染病的發(fā)生規(guī)律及預(yù)警分析中，SARIMA模型、ETS模型及ARIMA模型均是經(jīng)典、成熟的時間序列預(yù)測方法，也是我國疾病監(jiān)測工作者應(yīng)用最為廣泛的預(yù)測模型[12-13]。而HFRS作為具有典型季節(jié)性發(fā)病特征的傳染病，適用于SARIMA模型，能夠彌補既往研究中未納入HFRS季節(jié)性、周期性發(fā)病規(guī)律的ARIMA預(yù)測模型[4，13，16]的不足。ETS模型則利用序列以往實際觀測值的加權(quán)平均來預(yù)測未來數(shù)值，兼容全期平均和移動平均的優(yōu)勢，從遠到近給予呈指數(shù)形式逐漸減弱的權(quán)重[14]，同樣適用于預(yù)測HFRS的發(fā)病趨勢。本研究分別用SARIMA模型和ETS模型構(gòu)建湖南省HFRS病例預(yù)測模型，從模型的預(yù)測結(jié)果來看，SARIMA(1，0，0)(3，0，0)12預(yù)測結(jié)果平均絕對百分比誤差為13.18%，ETS模型預(yù)測結(jié)果平均絕對百分比誤差為31.14%，一般認為MAPE低于10%～15%時預(yù)測精度較好，SARIMA模型的預(yù)測精度要高于ETS模型。原因可能是ETS預(yù)測的前提是，歷史數(shù)據(jù)存在的各種因素的影響趨勢將持續(xù)保持，但是影響腎綜合征出血熱發(fā)病的因素是眾多的，當HFRS的發(fā)病情況出現(xiàn)上升或下降的趨勢時，指數(shù)平滑法就難以適應(yīng)。因此，SARIMA(1，0，0)(3，0，0)12模型更適合用于湖南省HFRS病例發(fā)病的短期預(yù)測[15]。

本研究所采用的SARIMA模型預(yù)測方法雖然較好地模擬和預(yù)測了湖南省HFRS發(fā)病趨勢，但無法揭示引起湖南省HFRS發(fā)病變化的各種影響因素[16]，同時所建模型是以歷史監(jiān)測數(shù)據(jù)序列為依據(jù)而建立的，僅適用于短期預(yù)測，仍需不斷用新的監(jiān)測數(shù)據(jù)對所建模型進行修正，并利用多種方法對影響HFRS發(fā)病的各種因素展開進一步研究[17]。