基于ARIMA乘積季節(jié)模型和Holt-Winters季節(jié)模型的梅毒月發(fā)病率預測

馬曉梅，史魯斌，其木格，閆國立，施學忠，孫春陽，徐學琴，趙倩倩

1)河南中醫(yī)藥大學公共衛(wèi)生與預防學科 鄭州 450046 2)河南省疾病預防控制中心免疫規(guī)劃科 鄭州 450046 3)鄭州大學公共衛(wèi)生學院衛(wèi)生統(tǒng)計學教研室 鄭州 450001

梅毒(syphilis)是由梅毒螺旋體引起的慢性、系統(tǒng)性性傳播疾病[1]。梅毒在全世界流行，近年來在我國的發(fā)病率有迅速升高的趨勢[1-2]。因此掌握其發(fā)展變化規(guī)律，早期預測其未來走勢，對有計劃地干預梅毒流行非常重要[1]。時間序列分析中的ARIMA乘積季節(jié)模型和Holt-Winters季節(jié)模型是預測某對象未來走勢的常用方法，技術(shù)成熟，短期預測精度高，在解釋預測變動規(guī)律方面具有重要價值，目前已廣泛應用于醫(yī)學衛(wèi)生領(lǐng)域[3-7]。大多數(shù)研究者[8-13]傾向于選擇ARIMA模型，因其側(cè)重綜合提取序列有效信息并能將其量化為具體函數(shù)表達式，但步驟復雜，結(jié)果難以解釋；Holt-Winters季節(jié)模型側(cè)重提取序列的確定性信息，目前在梅毒發(fā)病率預測中的應用價值尚待驗證[3,5-6]。作者以2005年1月至2016年6月梅毒月發(fā)病率數(shù)據(jù)為基礎(chǔ)，分別基于ARIMA乘積季節(jié)模型和Holt-Winters季節(jié)模型建模，探討最佳的梅毒發(fā)病預測方法，為梅毒防控提供依據(jù)。

1 資料與方法

1.1資料來源2005年1月至2016年6月全國梅毒月發(fā)病資料來自“疾病監(jiān)測信息報告管理系統(tǒng)”，全國人口資料來自《中國統(tǒng)計年鑒》。月發(fā)病率=梅毒月發(fā)病人數(shù)/人口數(shù)×105。

1.2建模方法以2005年1月至2015年12月梅毒月發(fā)病率數(shù)據(jù)為基礎(chǔ)，運用SPSS 22.0和Eviews 8.0軟件分別建立ARIMA乘積季節(jié)模型和Holt-Winters季節(jié)模型，采用2016年1至6月的實際數(shù)據(jù)驗證。選擇預測精度較高模型同法預測2016年7至12月梅毒月發(fā)病率。

1.3ARIMA乘積季節(jié)模型

1.3.1 模型介紹 ARIMA乘積季節(jié)模型稱為求和自回歸移動平均乘積季節(jié)模型，簡記為：ARIMA(p，d，q)×(P，D，Q)s，式中，p和q為自回歸和移動平均階數(shù)，d為差分次數(shù)，P和Q為季節(jié)性自回歸和移動平均階數(shù)，D為季節(jié)性差分次數(shù)，s為季節(jié)周期和循環(huán)長度[8]。

1.3.2 建模步驟 ①平穩(wěn)性檢驗。當序列非平穩(wěn)時，選擇合適的差分方式提取原序列中周期性和趨勢變動等確定性信息，確定d和D。②模型定階。利用樣本自相關(guān)系數(shù)與偏自相關(guān)系數(shù)的性質(zhì)確定p和q，而參數(shù)P和Q超過2階的情況很少見，可從低階到高階逐個嘗試[9]。③參數(shù)估計與模型檢驗。采用Box-Ljung統(tǒng)計量檢驗殘差序列是否為白噪聲序列。如果備選模型沒有通過，轉(zhuǎn)向步驟②，重新擬合。由于定階的主觀性導致有效模型并不唯一，故需在所有通過檢驗的模型中選擇BIC函數(shù)最小、擬合優(yōu)度較大、簡潔的模型為最優(yōu)模型。④序列預測。利用最優(yōu)模型預測，并用實際數(shù)據(jù)驗證。評價指標為預測誤差和平均絕對誤差(MAE)。其中，預測誤差=|預測值-實際值|/實際值×100%,MAE是所有單個預測誤差絕對值的平均值。

1.4Holt-Winters季節(jié)模型

1.4.1 模型介紹 Holt-Winters季節(jié)模型是一種指數(shù)平滑方法，可以同時修正序列的季節(jié)性和傾向性，包含Holt-Winters季節(jié)加法模型和Holt-Winters季節(jié)乘法模型。其基本思想是對同時具有趨勢效應和季節(jié)效應的序列分別進行指數(shù)平滑，目的是分解確定性信息，然后結(jié)合每種形式的平滑結(jié)果，對原時間序列進行預測。

1.4.2 模型方程 經(jīng)Holt-Winters季節(jié)乘法模型預測的時間序列yt平滑后的序列y't的方程式如下：

y't+k=at+btk+ct+k

a、b、c的計算公式如下：

at=α(yt-ct-s)+(1-a)(at-1+bt-1)

bt=β(at-at-1)+(1-β)bt-1

ct=γ(yt-at)+(1-γ)ct-s

at=α(yt/ct-s)+(1-a)(at-1+bt-1)

bt=β(at-at-1)+(1-β)bt-1

ct=γ(yt/at)+(1-γ)ct-s

上述式中，a表示整體平滑，b表示趨勢平滑，c為季節(jié)平滑;t=1,2，…,T；k為向后平滑期數(shù),k>0；α、β、γ為平滑參數(shù)，∈[0,1]，s為季節(jié)長度，對于月度數(shù)據(jù)s=12。

1.5統(tǒng)計工具采用Excel 2010 建立梅毒月發(fā)病率數(shù)據(jù)庫，運用SPSS 22.0和Eviews 8.0建立模型和分析數(shù)據(jù)，檢驗水準α=0.05。

2 結(jié)果

2.1梅毒月發(fā)病率特征分析見圖1。2005年1月至2016年6月我國梅毒月發(fā)病率(xt)有以下特點：①序列呈總體上升趨勢。②方差隨時間呈現(xiàn)較大波動，存在遞增型異方差。③序列蘊含以年為單位的周期性。我國梅毒月發(fā)病率呈現(xiàn)明顯季節(jié)性特征，發(fā)病高峰集中在7至8月，發(fā)病低谷集中在1至2月。

圖1 2005年1月至2016年6月梅毒月發(fā)病率時序圖

2.2ARIMA乘積季節(jié)模型建模

2.2.1 平穩(wěn)性檢驗 由于原序列前后方差波動較大，為消除異方差，先進行自然對數(shù)轉(zhuǎn)換，再進行1階12步差分提取趨勢和季節(jié)效應。變換后新序列(ΔΔ12lnxt)滿足建模平穩(wěn)性要求(圖2)。

圖2 ΔΔ12lnx序列的時序圖

2.2.2 模型定階 經(jīng)平穩(wěn)性檢驗后，可以確定s=12，d=1，D=1。ACF和PACF圖(圖3)顯示，自相關(guān)系數(shù)延遲1階后截尾，偏自相關(guān)系數(shù)拖尾，初步確定p=2，q=1；自相關(guān)系數(shù)延遲12階大于2倍標準差范圍，但是24階落入2倍標準差范圍，偏自相關(guān)系數(shù)顯示延遲12階、24階均在2倍標準差范圍，故P、Q嘗試0或1。

2.2.3 參數(shù)估計與模型檢驗 經(jīng)時間序列建模器反復調(diào)試、擬合與檢驗后，僅ARIMA(1,1,1)×(0,1,1)12和ARIMA(1,1,1)×(1,1,0)12參數(shù)全部有統(tǒng)計學意義(表1)。

表1 梅毒月發(fā)病率備選模型參數(shù)估計結(jié)果

對備選模型進行Box-Ljung檢驗，結(jié)果顯示兩者均通過白噪聲檢驗(P>0.05)。通過比較，ARI-MA(1,1,1)×(0,1,1)12模型BIC函數(shù)值較小，擬合優(yōu)度較大，模型簡潔，為相對最優(yōu)模型。綜上，確定模型口徑為：(1-B)(1-B12)(1+0.374B)xt=(1+0.740B)(1+0.775B12)εt，其中，B為延遲算子，εt為白噪聲序列。

2.2.4 序列預測 采用2016年1至6月梅毒實際月發(fā)病率驗證該模型預測結(jié)果(表2)。經(jīng)計算，MAE=2.704%。

2.3Holt-Winters季節(jié)模型建模及預測利用Eviews 8.0中的Exponential Smoothing，在平滑方法中分別選擇Holt-Winters Additive和Holt-Winters Multiplicative，平滑參數(shù)α、β和γ采用系統(tǒng)自動給定的方式產(chǎn)生。該模型預測值及預測誤差見表2。經(jīng)計算，Holt-Winters季節(jié)加法模型MAE=4.627%，Holt-Winters季節(jié)乘法模型MAE=3.794%。

表2 2016年1至6月梅毒月發(fā)病率預測值及預測誤差

2.4模型外推預測通過比較MAE，說明ARIMA乘積季節(jié)模型預測精度優(yōu)于Holt-Winters季節(jié)模型，故選用ARIMA乘積季節(jié)模型外推預測2016年7至12月梅毒月發(fā)病率、預測值及其95%CI,見表3。2005年至2016年預測值與實際值的比較見圖4。

表3 2016年7至12月梅毒月發(fā)病率ARIMA乘積季節(jié)模型預測值及95%CI

圖4 2005年至2016年梅毒月發(fā)病率實際值和預測值比較

3 討論

近年來，時間序列分析方法因其能根據(jù)歷史數(shù)據(jù)延伸預測而迅速活躍于醫(yī)學和公共衛(wèi)生領(lǐng)域[10-14]。該法以動態(tài)時間數(shù)據(jù)為基礎(chǔ)，通過分析時間序列的發(fā)展過程、方向和趨勢進行外推延伸，借以預測事物下一階段的水平。時間序列變動的規(guī)律性與不規(guī)律性共存，每一時期觀察值的大小均是各種因素在同一時刻發(fā)生作用的綜合結(jié)果，這些因素大致分為長期趨勢、周期性、隨機波動及綜合作用[3]。利用時間序列求出這些因素的數(shù)學模型后，可預測序列的未來走勢。當一個序列同時存在以上幾個因素時，提取最為常用的方法是ARIMA乘積季節(jié)模型和Holt-Winters季節(jié)模型。其中，ARIMA乘積季節(jié)模型是綜合性提取序列的發(fā)展特征，尤其針對長期趨勢、季節(jié)效應和隨機波動之間存在復雜交互關(guān)系的序列[3]。而Holt-Winters季節(jié)模型則是通過指數(shù)平滑的方法分解序列的變化規(guī)律，尤其針對確定性因素占明顯優(yōu)勢而不確定信息微弱的序列。其中，Holt-Winters季節(jié)加法模型適用于有線性趨勢且不依賴于序列水平的季節(jié)性效應的序列，Holt-Winters季節(jié)乘法模型適用于有線性趨勢且依賴于序列水平的季節(jié)性效應的序列[15]。

該研究利用2005年1月至2015年12月梅毒月發(fā)病率資料進行時間序列分析，分別建立ARIMA乘積季節(jié)模型和Holt-Winters季節(jié)模型。預測結(jié)果表明，ARIMA乘積季節(jié)模型預測誤差均在6%以下，MAE為2.704%；Holt-Winters季節(jié)加法模型有1組預測誤差為11.116%，其余在6%以下，MAE為4.627%；Holt-Winters季節(jié)乘法模型預測誤差均在6%以下，MAE為3.794%?？傮w來說，ARIMA乘積季節(jié)模型預測準確性較Holt-Winters季節(jié)模型高。其原因可能是ARIMA乘積季節(jié)模型不但能綜合且充分提取序列的長期趨勢、周期性和隨機波動等信息，而且能判斷這些因素之間確切的作用關(guān)系，其模型口徑為：(1-B)(1-B12)(1+0.374B)xt=(1+0.740B)(1+0.775B12)εt；而Holt-Winters季節(jié)模型僅能提取序列的確定性信息，對隨機波動信息浪費嚴重。但從建模實際操作過程來看，ARIMA乘積季節(jié)模型建模過程比較繁瑣，在模型定階的選擇上存在多種評判標準，結(jié)果難以解釋。Holt-Winters季節(jié)模型操作較為簡便，3個重要的平滑參數(shù)α、β和γ均經(jīng)系統(tǒng)多次比較，自動產(chǎn)出最優(yōu)參數(shù)，結(jié)果較易解釋。故當序列中確定性信息較強時，選擇Holt-Winters季節(jié)模型預測也會得到非常不錯的預測結(jié)果[3]。

ARIMA乘積季節(jié)模型和Holt-Winters季節(jié)模型均屬于短期預測精度較高的模型。隨著時間的推移，未知信息增多，模型估計精度也會逐漸變差[3]。該研究所采用的數(shù)據(jù)是全國梅毒月發(fā)病率數(shù)據(jù)，在預測時僅僅是從序列本身出發(fā)，沒有考慮經(jīng)濟、社會等影響因素，因此可能會出現(xiàn)一定的預測誤差。故在實際操作中，應不斷更新數(shù)據(jù)對模型進行重新擬合，以提高預測的準確性和適用性，從而有效指導我國的梅毒防控工作。

[1] ZHANG X,ZHANG T,PEI J,et al.Time series modelling of syphilis incidence in China from 2005 to 2012[J].PLoS One,2016,11(2):e0149401

[2] 王小麗,楊永利,施學忠,等.幾種預測模型對中國梅毒發(fā)病率預測結(jié)果的比較[J].鄭州大學學報(醫(yī)學版),2015,50(2):164

[3] 王燕.應用實踐序列分析[M].3版.北京:中國人民大學出版社,2014:160

[4] LIU L,LUAN RS,YIN F,et al.Predicting the incidence of hand, foot and mouth disease in Sichuan province, China using the ARIMA model[J].Epidemiol Infect,2016,144(1):144

[5] 陳媛.應用holt-winters加法模型預測出院人次[J].中國衛(wèi)生統(tǒng)計,2012,29(2):260

[6] 史蕓萍，馬家奇.指數(shù)平滑法在流行性腮腺炎預測預警中的應用[J].中國疫苗和免疫，2010,16(3)：233

[7] CHAYA S,DANGOR Z,SOLOMON F,et al.Incidence of tuberculosis meningitis in a high HIV prevalence setting: time-series analysis from 2006 to 2011[J].Int J Tuberc Lung Dis,2016,20(11):1457

[8] 楊小兵,孔德廣,江高峰.ARIMA乘積季節(jié)模型在手足口病發(fā)病預測中的應用研究[J].中國預防醫(yī)學雜志,2016,17(3):207

[9] 于林鳳,吳靜,周鎖蘭,等.ARIMA季節(jié)模型在我國丙肝發(fā)病預測中的應用[J].鄭州大學學報(醫(yī)學版),2014,49(3):344

[10]SONG X,XIAO J,DENG J,et al.Time series analysis of influenza incidence in Chinese provinces from 2004 to 2011[J].Medicine(Baltimore),2016,95(26):e3929

[11]CHEN BH,SUMI A,TOYODA S,et al.Time series analysis of reported cases of hand, foot, and mouth disease from 2010 to 2013 in Wuhan, China[J].BMC Infect Dis,2015,15(1):1

[12]WANG T,ZHOU Y,WANG L,et al.Using an autoregressive integrated moving average model to predict the incidence of hemorrhagic fever with renal syndrome in Zibo, China, 2004-2014[J].Jpn J Infect Dis,2016,69(4):279

[13]萬燕麗,楊永利,施念,等.ARIMA模型在河南省 AIDS疫情預測中的應用[J].鄭州大學學報(醫(yī)學版),2015,50(2):160

[14]胡建利,劉文東,梁祁,等.季節(jié)指數(shù)法和ARIMA模型在感染性腹瀉周發(fā)病數(shù)預測中的應用研究[J].中華疾病控制雜志,2013,17(8):718

[15]武松.SPSS統(tǒng)計分析大全[M].北京:清華大學出版社,2014:358