祝寒松,黃文龍,謝忠杭,陳光敏,曹 洋,章?tīng)N明,陳 武,歐劍鳴,洪榮濤
SARIMA模型在戊型肝炎發(fā)病趨勢(shì)預(yù)警預(yù)測(cè)中的應(yīng)用
祝寒松1,黃文龍1,謝忠杭1,陳光敏1,曹 洋2,章?tīng)N明1,陳 武1,歐劍鳴1,洪榮濤1
目的 采用SARIMA模型對(duì)福建省戊肝發(fā)病趨勢(shì)進(jìn)行預(yù)測(cè),為預(yù)警和風(fēng)險(xiǎn)評(píng)估提供定量數(shù)據(jù)。方法 以Eviews5.0對(duì)福建省2004年1月-2015年12月戊肝的月發(fā)病例數(shù)進(jìn)行SARIMA分析。結(jié)果 2004年1月-2014年8月福建省戊肝月發(fā)病序列呈先升后降的趨勢(shì)和周期性波動(dòng),取自然對(duì)數(shù)和1次1階非季節(jié)差分后序列得到平穩(wěn)化,模型SARIMA (0,1,1)(1,0,1)12和SARIMA(2,1,0)(1,0,1)12參數(shù)有統(tǒng)計(jì)學(xué)意義,殘差為白噪聲,后者為最優(yōu)模型,表達(dá)式為:(1+0.61L+0.23L2)(1-0.89L12)(1-L)log(yt)=(1-0.89L12)εt。靜態(tài)回代預(yù)測(cè)值和實(shí)際值吻合較好,相對(duì)誤差取絕對(duì)值后均數(shù)為13.06%。2014年1-8月預(yù)測(cè)值的相對(duì)誤差較小(5月份除外),預(yù)測(cè)標(biāo)準(zhǔn)誤(S.E)較小。結(jié)論 運(yùn)用SARIMA模型可對(duì)戊肝發(fā)病趨勢(shì)進(jìn)行較準(zhǔn)確的短期預(yù)測(cè),可為及時(shí)、科學(xué)地研判風(fēng)險(xiǎn)提供可靠的數(shù)據(jù)基礎(chǔ)。
季節(jié)時(shí)間序列模型;戊肝;預(yù)警;預(yù)測(cè);風(fēng)險(xiǎn)評(píng)估
戊型肝炎(戊肝,Hepatitiv E)是由戊型肝炎病毒引起的一種急性人獸共患病,目前在全球各地廣泛傳播,嚴(yán)重危害人類(lèi)健康,據(jù)估計(jì)20億的世界人口有可能被感染Hepatitiv E病毒的風(fēng)險(xiǎn),死亡率也占0.5%,并且有可能在許多發(fā)展中國(guó)家導(dǎo)致大暴發(fā),近年來(lái)發(fā)病率在我國(guó)呈逐年上升趨勢(shì)[1-3]。為了對(duì)戊肝發(fā)病趨勢(shì)進(jìn)行預(yù)測(cè)和預(yù)警,考慮到其發(fā)病具有一定的周期性,而季節(jié)時(shí)間序列模型(Seasonal Autoregressive Integrated Moving Average Model, SARIMA)可用于對(duì)存在周期性變化的時(shí)間序列進(jìn)行預(yù)測(cè),因此,本文擬采用SARIMA模型對(duì)福建省戊肝發(fā)病趨勢(shì)進(jìn)行定量預(yù)測(cè),以期為風(fēng)險(xiǎn)評(píng)估提供基礎(chǔ)資料,現(xiàn)報(bào)告如下。
1.1 材料 來(lái)自福建省各級(jí)衛(wèi)生機(jī)構(gòu)通過(guò)中國(guó)疾病預(yù)防控制信息系統(tǒng)的報(bào)告(按“現(xiàn)住址”、“終審”選項(xiàng)下載)。
1.2 方法 采用SAS9.2軟件對(duì)數(shù)據(jù)進(jìn)行匯總、清洗,運(yùn)用Eviews 5.0軟件對(duì)數(shù)據(jù)進(jìn)行SARIMA建模和分析,α=0.05(雙側(cè))。
1.2.1 原理 ARIMA模型通常借助時(shí)間序列的隨機(jī)特性來(lái)描述事物的發(fā)展變化規(guī)律,從而解釋并預(yù)測(cè)時(shí)間序列的變化發(fā)展規(guī)律,而SARIMA模型是ARIMA模型的擴(kuò)展和改進(jìn),包含了季節(jié)性或周期性因素[4-5],對(duì)于存在季節(jié)性的非平穩(wěn)時(shí)間序列不能直接建立ARIMA模型,可考慮進(jìn)行非季節(jié)差分和季節(jié)性差分使其平穩(wěn)化,基本模型:SARIMA(p,d,q)(P,D,Q)S,表達(dá)式為:φp(L)ΦP(LS)(1-L)d(1-LS)Dyt=θq(L)ΘQ(LS)εt。其中,p、q分別是非季節(jié)自回歸過(guò)程AR和移動(dòng)平均過(guò)程MA的階數(shù);d、D分別是序列yt的非季節(jié)差分和季節(jié)差分階數(shù);P、Q分別是季節(jié)自回歸過(guò)程SAR和移動(dòng)平均過(guò)程SMA的階數(shù)。φp(L)、ΦP(L)分別是非季節(jié)自回歸過(guò)程AR和季節(jié)自回歸過(guò)程SAR的滯后算子多項(xiàng)式,θq(L)、ΘQ(LS)分別是非季節(jié)移動(dòng)平均過(guò)程MA和季節(jié)移動(dòng)平均過(guò)程SMA的滯后算子多項(xiàng)式,(1-L)d、(1-LS)D分別是對(duì)序列yt的非季節(jié)差分和季節(jié)差分滯后算子,S是季節(jié)差分的步長(zhǎng),εt是殘差序列(白噪聲序列)。
1.2.2方法 ①繪制原序列的曲線圖以識(shí)別其基本形式,加入截距項(xiàng)和時(shí)間趨勢(shì)項(xiàng)做ADF單位根檢驗(yàn)和Q統(tǒng)計(jì)量以判斷平穩(wěn)性。若為非平穩(wěn),可考慮對(duì)其取自然對(duì)數(shù)后做非季節(jié)差分和(或)季節(jié)性差分使之滿足平穩(wěn)性條件。②識(shí)別新序列自相關(guān)函數(shù)(autocorrelation function, ACF)和偏相關(guān)函數(shù)(partial autocorrelation function, PACF) 的SARIMA模型形式。③估計(jì)模型參數(shù),對(duì)結(jié)果進(jìn)行t檢驗(yàn)(參數(shù)檢驗(yàn)),通過(guò)Q檢驗(yàn)進(jìn)行殘差分析(白噪聲檢驗(yàn))。如果殘差序列不是白噪聲序列,說(shuō)明殘差序列中還存在有用的信息未被提取出來(lái),需要對(duì)原模型進(jìn)一步調(diào)整,以便得到更合適的模型。④利用調(diào)整R2、Akaike info criterion (AIC)準(zhǔn)則和Schwarz criterion (SC)準(zhǔn)則評(píng)價(jià)其相對(duì)優(yōu)勢(shì),調(diào)整R2越大、AIC和SC越小,一般認(rèn)為越好。⑤利用所估計(jì)的SARIMA模型,進(jìn)行回代評(píng)價(jià)和預(yù)測(cè)。
2.1 原序列分析 2004年1月-2014年8月福建省戊肝月發(fā)病序列呈先升后降的趨勢(shì)和周期性波動(dòng)(圖1),ADF單位根檢驗(yàn),P=0.00<0.05,但自相關(guān)函數(shù)呈現(xiàn)指數(shù)衰減,且衰減速度緩慢,所以認(rèn)為該序列不平穩(wěn)(圖2)。
圖1 福建省2004-2015年戊肝月發(fā)病例數(shù)的實(shí)際值和預(yù)測(cè)值
Fig.1 Observed and predicted values of hepatitis E in Fujian Province, per month from 2004 to 2015
圖2 福建省2004年1月-2014年8月戊肝月發(fā)病數(shù)序列的自相關(guān)和偏相關(guān)圖
Fig.2 ACF and PACF diagram of hepatitis E incidence in Fujian Province, per month from January 2004 to August 2014
2.2 平穩(wěn)化 采用2種方法對(duì)原序列進(jìn)行處理,第1種是取自然對(duì)數(shù)和1次1階非季節(jié)差分,再做1次12階季節(jié)差分;第2種是只取自然對(duì)數(shù)和1次1階非季節(jié)差分。通過(guò)比較,后者效果更好,現(xiàn)只對(duì)其進(jìn)行分析:差分后的序列已接近平穩(wěn)(圖3),ADF單位根檢驗(yàn)P=0.00<0.05;ACF和PACF顯示,序列平穩(wěn)性已得到明顯改進(jìn),滯后1階和7階顯著地不為0,其它近似為1個(gè)平穩(wěn)過(guò)程(圖4)。因此,可認(rèn)為序列已得到平穩(wěn)化。
圖3 經(jīng)自然對(duì)數(shù)和1次1階非季節(jié)差分后的序列
Fig.3 Time series diagram after conducting natural logarithm and once non-seasonal of lag 1 difference
圖4 序列經(jīng)差分后的自相關(guān)和偏相關(guān)圖
Fig.4 ACF and PACF diagram after conducting natural logarithm and a non-seasonal difference
2.3 建模 由2.2差分可知d=1、D=0,ACF滯后1階呈截尾,PACF滯后2階呈截尾或者拖尾至2階。通過(guò)嘗試,2組模型通過(guò)t檢驗(yàn)(p<0.05)和Q檢驗(yàn)(全部p>0.05,已呈明顯的白噪聲序列):SARIMA(0,1,1)(1,0,1)12、SARIMA(2,1,0) (1,0,1)12(表1)。通過(guò)調(diào)整R2、AIC和SC比較,SARIMA(2,1,0) (1,0,1)12擬合優(yōu)度較好(表2),Q檢驗(yàn)見(jiàn)圖5,表達(dá)式為:(1+0.61L+0.23L2)(1-0.89L12)(1-L)log(yt)=(1-0.89L12)εt。
圖5 SARIMA(2,1,0)(1,0,1)12模型的Q檢驗(yàn)
變量VariableSARIMA(0,1,1)(1,0,1)12系數(shù)CoefficienttPSARIMA(2,1,0)(1,0,1)12系數(shù)CoefficienttPar(1)----0.61-6.470.00ar(2)----0.23-2.480.01sar(12)0.8633.280.000.8939.450.00ma(1)-0.40-5.480.00---sma(12)-0.87-28.120.00-0.90-39.070.00
表2 2組SARIMA模型的擬合優(yōu)度比較
2.4 預(yù)測(cè) 采用SARIMA(2,1,0)(1,0,1)12模型回代靜態(tài)預(yù)測(cè)評(píng)估,擬合結(jié)果顯示,預(yù)測(cè)值和實(shí)際值吻合較好,2004年1月-2005年3月的實(shí)際值經(jīng)差分已去除(圖1),2005年4月-2014年8月相對(duì)誤差均數(shù)為0.34%,取絕對(duì)值后相對(duì)誤差均數(shù)為13.06%。2014年1-8月預(yù)測(cè)值的相對(duì)誤差較小(5月份除外),預(yù)測(cè)標(biāo)準(zhǔn)誤(S.E)較小(表3)。2014年9月-2015年12月進(jìn)行動(dòng)態(tài)預(yù)測(cè),S.E逐漸較大(圖1),2014年9-12月的S.E和95%CI上下限見(jiàn)表3。
表3 福建省2014年戊肝預(yù)測(cè)情況
戊肝傳播途徑多樣,其中以糞-口途徑傳播為主,而且有證據(jù)顯示通過(guò)攝食受感染的肉類(lèi)可以引起人獸共患傳播,易引起食物或水源性流行暴發(fā),另外,普通人群戊肝病毒的感染率較高,提示亞臨床感染的存在,也是無(wú)癥狀感染、重癥肝炎、不明原因肝炎以及暴發(fā)性急性肝炎的重要原因之一[6-7],然而,目前該病尚缺乏特效的治療方法,也沒(méi)有特異性免疫制劑可供預(yù)防,而福建省省戊肝發(fā)病總體呈上升態(tài)勢(shì)。因此,加強(qiáng)戊肝監(jiān)測(cè)資料分析并予以預(yù)測(cè)預(yù)警以嚴(yán)防流行和暴發(fā)顯得尤為重要。
但傳染病的暴發(fā)流行受到多種不確定因素的影響,這使得在傳染病早期的預(yù)測(cè)預(yù)警上存在著諸多困難,導(dǎo)致了傳染病早期預(yù)防控制工作一直較為滯后,因此如何及時(shí)有效地預(yù)測(cè)預(yù)警傳染病的暴發(fā)流行一直是傳染病預(yù)測(cè)控制工作的重點(diǎn)[8-10]。
應(yīng)用模型對(duì)傳染病的發(fā)病及流行強(qiáng)度進(jìn)行預(yù)測(cè)預(yù)警有多種方法:灰色動(dòng)態(tài)模型GM(1,1)將原始序列累加、求均值而生成新的數(shù)列,使得GM(1,1)模型預(yù)測(cè)精度降低;利用人工神經(jīng)網(wǎng)絡(luò)(ANN)模型預(yù)測(cè)發(fā)病率,其難度在于確定網(wǎng)絡(luò)結(jié)構(gòu),隱含層節(jié)點(diǎn)數(shù)太少,預(yù)測(cè)精度無(wú)法保證,節(jié)點(diǎn)數(shù)太多,又易陷入局部極小值,因此如何選擇一個(gè)最佳的網(wǎng)絡(luò)結(jié)構(gòu),成為一個(gè)關(guān)鍵問(wèn)題;傳統(tǒng)的時(shí)間序列模型要求序列具有平穩(wěn)的線性趨勢(shì),但許多傳染病的發(fā)病情況有著明顯的季節(jié)性和周期性,如果不考慮這些因素的影響,做出的預(yù)測(cè)往往不精確[11-12]。而SARIMA模型適用于多種復(fù)雜的時(shí)間序列模式,可將擬合誤差作為重要因素納入模型中,因此預(yù)測(cè)精度較高,作為疾病風(fēng)險(xiǎn)評(píng)估、預(yù)防和控制的監(jiān)測(cè)數(shù)據(jù),這對(duì)于規(guī)劃公共衛(wèi)生干預(yù)措施是一個(gè)非常有價(jià)值的工具[13-14]。
福建省戊肝月發(fā)病序列總體呈現(xiàn)明顯的周期性波動(dòng)(圖1),而且自相關(guān)函數(shù)呈指數(shù)緩慢衰減,為非平穩(wěn)性,因此該序列比較適合建立SARIMA模型定量預(yù)測(cè)月(周)發(fā)病數(shù)(率),可以為傳染病風(fēng)險(xiǎn)評(píng)估提供基礎(chǔ)數(shù)據(jù)。
分析顯示,SARIMA(2,1,0)(1,0,1)12模型擬合效果較好,回代預(yù)測(cè)值和實(shí)際值相當(dāng)吻合(平均誤差為13.06%),對(duì)2014年的預(yù)測(cè)標(biāo)準(zhǔn)誤較低,表明該模型預(yù)測(cè)較為準(zhǔn)確、精度較高,可以認(rèn)為它對(duì)于原序列是個(gè)理想的模型。但是也要注意到,圖1中紫色和綠色間的2014-2015年預(yù)測(cè)置信區(qū)間逐漸變寬,從而表明預(yù)測(cè)期越往后,模型的預(yù)測(cè)精度越差,因此比較適用于短期預(yù)測(cè)。在實(shí)際應(yīng)用中,應(yīng)該不斷加入新的實(shí)際值,然后修正模型和重新擬合預(yù)測(cè)值[15],所以對(duì)于2015年,在2014年的實(shí)際值納入序列矯正之后再做以預(yù)測(cè),結(jié)果才能更可靠。
值得一提的是,除5月份外,2014年1-8月的相對(duì)誤差較小,是較為成功的模型預(yù)測(cè)。至于5月份,其實(shí)際值偏低,且低于預(yù)測(cè)值的95%CI下限,同時(shí)大幅低于2011-2013年同期實(shí)際值(113、104和124例),需要從網(wǎng)絡(luò)直報(bào)、自然環(huán)境因素和防控力度的加強(qiáng)等多方面因素考慮。另外,本研究顯示,2014年2-4月和7、8月實(shí)際值均不高于95%CI上限預(yù)警值,實(shí)際情況也無(wú)流行和暴發(fā),與預(yù)測(cè)基本一致。
但也應(yīng)注意到,由于SARIMA模型也有不足之處,一是建立合理的模型需要大量的觀測(cè)數(shù)據(jù)(至少50個(gè)樣本);二是該模型是基于正態(tài)假設(shè),對(duì)于罕見(jiàn)病,即便是對(duì)數(shù)或平方根轉(zhuǎn)換后也不是十分合適。所以對(duì)于小樣本,可以結(jié)合泊松回歸模型預(yù)測(cè),但對(duì)于大樣本,前者的預(yù)測(cè)精確度高于后者[16]。
[1]Lin J, Norder H, Uhlhorn H, et al. Novel hepatitis E like virus found in Swedish moose[J]. J Gen Virol, 2014, 95(Pt 3): 557-570. DOI: 10.1099/vir.0.059238-0
[2]Kuniholm MH, Purcell RH, Mc Quillan GM, et al. Epidemiology of hepatitis E virus in the United States: results from the Third National Health and Nutrition Examination Survey, 1988-1994[J]. J Infect Dis, 2009, 200(1): 48-56. DOI: 10.1086/599319
[3]Wei S, Lu YH, Gao MY, et al. Time series analysis of hepatitis E incidence in China[J]. Chin J Health Statistics, 2012, 29(6): 808-811. (in Chinese) 魏珊, 陸一涵, 高眉揚(yáng), 等. 我國(guó)戊型肝炎發(fā)病例數(shù)的時(shí)間序列分析[J]. 中國(guó)衛(wèi)生統(tǒng)計(jì), 2012, 29(6): 808-811.
[4]Moosazadeh M, Nasehi M, Bahrampour A, et al. Forecasting tuberculosis incidence in Iran using Box-Jenkins models[J]. Iran Red Crescent Med J, 2014, 16(5): 1-6. DOI: 10.5812/ircmj.11779
[5]Kam HJ, Sung JO, Park RW. Prediction of daily patient numbers for a regional emergency medical center using time series analysis[J]. Healthc Inform Res, 2010, 16(3): 158-165. DOI: 10.4258/hir.2010.16.3.158
[6]Wang HR, Yan YS, Xiao JX, et al. Seroepidemiological investigation and analysis of hepatitis virus infection among various groups of population in Fujian province[J]. Chin J Zoonoses, 2007, 23(4): 370-372. (in Chinese) 王惠榕, 嚴(yán)延生, 蕭劍雄, 等. 福建省不同人群中戊型肝炎病毒感染的血清流行病學(xué)調(diào)查分析[J]. 中國(guó)人獸共患病學(xué)報(bào), 2007, 23(4): 370-372
[7]Arends JE, Ghisetti V, Irving W, et al. Hepatitis E: An emerging infection in high income countries[J].J Clin Virol, 2014, 59(2): 81-88. DOI: 10.1016/j.jcv.2013.11.013
[8]Yang Z, Ye ZH, You AG, et al. Application of multiple seasonal ARIMA model in prediction of tuberculosis incidence[J]. Chin J Public Health, 2013, 29(4): 469-473. (in Chinese) 楊召, 葉中輝, 尤愛(ài)國(guó), 等. 乘積季節(jié)ARIMA模型在結(jié)核病發(fā)病率預(yù)測(cè)中的應(yīng)用[J]. 中國(guó)公共衛(wèi)生, 2013, 29(4): 469-473.
[9]Song Q. On the weight convergence of Elman networks[J]. IEEE Trans Neural Netw, 2010, 21(3): 463-480. DOI: 10.1109/TNN.2009.2039226
[10]Lai SJ, Li ZJ, Jin LM, et al. Evaluation content and their indicators of early warning system for infectious disease outbreak[J]. Chin J Epidemiol, 2009, 30(6): 637-641. (in Chinese) 賴圣杰, 李中杰, 金連梅, 等. 傳染病暴發(fā)早期預(yù)警系統(tǒng)評(píng)價(jià)內(nèi)容及其指標(biāo)[J]. 中華流行病學(xué)雜志, 2009, 30(6): 637-641.
[11]Bras AL, Gomes D, Filipe PA, et al. Trends, seasonality and forecasts of pulmonary tuberculosis in Portugal[J]. Int J Tuberc Lung Dis, 2014, 18(10): 1202-1210. DOI: 10.5588/ijtld.14.0158
[12]Shen TQ, Liu WD, Hu JL, et al. The application of x-11-ARIMA process in dysentery prediction[J]. Chin J Health Statistics, 2014, 31(3): 395-398. (in Chinese) 申銅倩, 劉文東, 胡建立, 等. X-11-ARIMA過(guò)程在痢疾疫情預(yù)測(cè)中的應(yīng)用研究[J]. 中國(guó)衛(wèi)生統(tǒng)計(jì),2014, 31(3): 395-398.
[13]An SY, Zhao Z, Guo JQ, et al. Forecasting measles epidemic situation by applying the time series model in Liaoning Province[J]. Chin J Health Statistics, 2014, 31(5): 781-783. (in Chinese) 安淑一, 趙卓, 郭軍巧, 等. 應(yīng)用時(shí)間序列模型預(yù)測(cè)遼寧省麻疹疫情[J]. 中國(guó)衛(wèi)生統(tǒng)計(jì), 2014, 31(5): 781-783.
[14]Feng HF, Duan GC, Zhang RG, et al. Time series analysis of hand-foot-mouth disease hospitalization in Zhengzhou: establishment of forecasting models using climate variables as predictors[J]. PLoS One, 2014, 9(1): 1-10. DOI: 10.1371/journal.pone.0087916
[15]Mehdi K, Mehdi B, Heiazi SR. Combining seasonal ARIMA models with computational intelligence techniques for time series forecasting[J]. Soft Comput, 2012, 16(6): 1091-1105. DOI: 10.1007/s00500 -012-0805-9
[16]Hu WB, Tong SL, Mengersen K, et al. Weather variability and the incidence of cryptosporidiosis: comparison of time series poisson regression and SARIMA models[J]. Elsevier Inc, 2007, 17(9): 679-688. DOI: 10.1016/j.annepidem.2007.03.020
SARIMA model for incidence trend and prediction of hepatitis E
ZHU Han-song1,HUANG Wen-long1,XIE Zhong-hang,CHEN Guang-min,CAO Yang2,ZHANG Can-ming1,CHEN Wu1,OU Jian-ming1,HONG Rong-tao1
(1.FujianProvincialCenterforDiseaseControl&Prevention,Fuzhou350001,China;2.ChineseCenterforDiseaseControlandPrevention,Bejing102206,China)
We forecasted the incidence trends of hepatitis E in Fujian Province by SARIMA model, in order to provide quantitative data for early warning and risk assessment. The monthly cases of hepatitis E in Fujian Province from January 2004 to December 2015 were analyzed for SARIMA using Eviews 5.0 software. Results showed that monthly incidence sequences of hepatitis E in Fujian Province, from January 2004 to August 2014, showed a trend of increased and then decreased, as well as cyclical fluctuations. And it was into smoothed after done natural logarithm and once non-seasonal of lag 1 difference. SARIMA (0,1,1)(1,0,1)12and SARIMA(2,1,0) (1,0,1)12were statistically significant,and the residual was white noise, in which the latter was the optimal model, expressed as: (1+0.61L+0.23L2)(1-0.89L12)(1-L)log(yt)=(1- 0.89L12)εt. Static back generation forecast and actual values were in good agreement, in which the mean of relative error was 13.06% after taking the absolute value. Relative error of predictive value in 2014 January to August (except in May) was small, and the prediction standard error (S.E) small, too. Incidence trends of hepatitis E could short-term predicted through SARIMA model, thus it could provide reliable data base in order to judge the risk of infectious disease more timely and scientifically.
SARIMA; hepatitis E; warning; forecast; risk assessment
Horg Rong-tao, Email:hrt@fjcdc.com.cn
洪榮濤,Email:hrt@fjcdc.com.cn
1.福建省疾病預(yù)防控制中心,福建省人獸共患病研究重點(diǎn)實(shí)驗(yàn)室,福州 350001; 2.中國(guó)疾病預(yù)防控制中心,應(yīng)急中心監(jiān)測(cè)預(yù)警與風(fēng)險(xiǎn)評(píng)估辦公室,北京 102206
10.3969/cjz.j.issn.1002-2694.2015.02.014
R373
A
1002-2694(2015)02-0158-05
2014-09-12;
2014-11-23