于林鳳，吳 靜，周鎖蘭，丁 勇#

1)南京醫(yī)科大學(xué)生物醫(yī)學(xué)工程系 南京 210029 2)南京醫(yī)科大學(xué)數(shù)學(xué)與計(jì)算機(jī)教研室 南京 210029

ARIMA季節(jié)模型在我國丙肝發(fā)病預(yù)測中的應(yīng)用*

于林鳳1)，吳 靜2)，周鎖蘭1)，丁 勇2)#

1)南京醫(yī)科大學(xué)生物醫(yī)學(xué)工程系 南京 210029 2)南京醫(yī)科大學(xué)數(shù)學(xué)與計(jì)算機(jī)教研室 南京 210029

#通訊作者，男，1956年8月生，碩士，教授，研究方向：生物統(tǒng)計(jì)，E-mail:yding@njmu.edu.cn

ARIMA季節(jié)模型；丙肝；發(fā)??；預(yù)測

目的：應(yīng)用ARIMA季節(jié)模型對我國丙肝發(fā)病進(jìn)行預(yù)測。方法利用2004年至2011年我國丙肝的月發(fā)病數(shù)建立ARIMA季節(jié)模型，對2012年丙肝的月發(fā)病數(shù)進(jìn)行預(yù)測，并用實(shí)際數(shù)據(jù)評(píng)估模型的預(yù)測效果。同法對同期甲肝發(fā)病數(shù)據(jù)進(jìn)行建模和預(yù)測。對丙肝和甲肝2004年至2011年的月發(fā)病數(shù)按年歸一化處理后計(jì)算方差。比較甲肝和丙肝的預(yù)測效果。結(jié)果成功建立ARIMA(1,1,1)(2,1,0)12季節(jié)模型，模型的表達(dá)式為：(1+0.222L)(1+0.820L12+0.694L24)(1-L)(1-L12)lnYt=(1+0.648L)εt,參數(shù) AR(1)=-0.222(t=-2.392,P=0.020)，SAR(12)=-0.820(t=-8.009,P<0.001),SAR(24)=-0.694(t=-6.124,P<0.001),MA(1)=-0.648(t=-5.889,P<0.001),殘差序列是白噪聲序列(P>0.05)；模型擬合效果的R2為0.824，預(yù)測的平均相對誤差為0.078。歸一化后丙肝和甲肝發(fā)病數(shù)的平均方差分別為0.030和0.047，提示丙肝原始數(shù)據(jù)周期性動(dòng)態(tài)變化較甲肝更趨一致。甲肝預(yù)測的平均相對誤差為0.138，大于丙肝。結(jié)論ARIMA(1,1,1)(2,l,0)12季節(jié)模型可用于預(yù)測我國丙肝的發(fā)病規(guī)律。樣本數(shù)據(jù)的周期性動(dòng)態(tài)變化趨勢越一致，ARIMA季節(jié)模型的預(yù)測結(jié)果也越準(zhǔn)確。

丙型病毒性肝炎(簡稱為丙肝)是一種由丙型肝炎病毒感染引起的病毒性肝炎，是一種對人群健康及生命危害較大的疾病，50%～80%的感染者進(jìn)展為慢性狀態(tài)，可導(dǎo)致肝臟慢性炎癥壞死及纖維化，其中20%～30%的患者發(fā)展為肝硬化甚至肝細(xì)胞癌。丙肝是歐美及日本等國家終末期肝病的最主要原因之一。我國衛(wèi)生部《2011年度全國法定傳染病報(bào)告發(fā)病、死亡統(tǒng)計(jì)表》數(shù)據(jù)顯示，2011年我國報(bào)告的丙肝發(fā)病數(shù)量超過17萬例，比2010年增長了13.1%，是乙肝的5倍。我國多地已出現(xiàn)丙肝疫情的暴發(fā)，丙肝的防控形勢相當(dāng)嚴(yán)峻[1]。

探討疾病的流行規(guī)律、預(yù)測發(fā)病人數(shù)對傳染性疾病的防治工作有著重要的指導(dǎo)意義。求和自回歸移動(dòng)平均(autoregressive integrated moving average, ARIMA)模型是一種基于時(shí)間序列分析、預(yù)測和控制的方法，其基本思想是利用時(shí)間序列的觀測值所具有的依存關(guān)系或自相關(guān)性，預(yù)測對象發(fā)展的未來值。隨機(jī)季節(jié)模型與ARIMA模型的結(jié)合即為ARIMA季節(jié)模型[2-4]，該模型能綜合考慮季節(jié)、趨勢和隨機(jī)干擾等因素，預(yù)測效果較好。作者采用ARIMA季節(jié)模型對我國丙肝的發(fā)病進(jìn)行了預(yù)測，并探討了數(shù)據(jù)的周期性變化趨勢與模型預(yù)測效果的關(guān)系。

1 資料與方法

1.1資料來源丙肝和甲肝數(shù)據(jù)資料來源于我國衛(wèi)生部網(wǎng)站(http://www.moh.gov.cn)2004年1月至2012年12月的全國法定報(bào)告?zhèn)魅静∫咔橘Y料，其中2004年1月至2011年12月的數(shù)據(jù)用于建立模型，2012年1月至12月的數(shù)據(jù)用于驗(yàn)證模型的預(yù)測效果。

1.2建模方法ARIMA季節(jié)模型一般表示為ARIMA(p,d,q)(P,D,Q)s，其中各參數(shù)的意義為：非季節(jié)差分階數(shù)d和季節(jié)差分階數(shù)D，非季節(jié)自回歸階數(shù)p和季節(jié)自回歸階數(shù)P，非季節(jié)移動(dòng)平均階數(shù)q和季節(jié)移動(dòng)平均階數(shù)Q，季節(jié)的長度s。當(dāng)P=D=0時(shí)，該模型便是一般的ARIMA模型。此次研究采用Eviews 6.0進(jìn)行建模和數(shù)據(jù)的處理及分析[5-6]。建模過程如下。

①序列的平穩(wěn)化：在確定時(shí)間序列模型之前需把不平穩(wěn)的時(shí)間序列轉(zhuǎn)化為平穩(wěn)的序列。為消除異方差，首先對原始數(shù)據(jù)進(jìn)行自然對數(shù)轉(zhuǎn)換，以平穩(wěn)序列的方差，然后根據(jù)變換后序列的自相關(guān)(ACF)和偏自相關(guān)(PACF)圖，確定d和D，d和D宜取較低階。s根據(jù)疾病的背景知識(shí)獲得。選擇不含常數(shù)項(xiàng)和趨勢項(xiàng)的增廣的迪基福勒檢驗(yàn)法(augmented Dickey-Fuller test,ADF)對結(jié)果進(jìn)行分析，檢驗(yàn)水準(zhǔn)α=0.05。

②模型的識(shí)別：對于變換后的平穩(wěn)時(shí)間序列，觀察其ACF和PACF圖，確定p、q值。參數(shù)P、Q超過2階的情況很少見[3]，可以分別取0、1、2并由低階到高階逐個(gè)試驗(yàn)，根據(jù)模型的擬合優(yōu)度、殘差情況以及系數(shù)間的相關(guān)性進(jìn)行綜合判斷，以確定這兩個(gè)參數(shù)。

③參數(shù)估計(jì)及檢驗(yàn)：運(yùn)用最大似然法或無約束最小二乘法估計(jì)模型的系數(shù)，并進(jìn)行參數(shù)估計(jì)和檢驗(yàn)，用擬合優(yōu)度比較該模型與其他模型的優(yōu)劣。

④模型的診斷：一個(gè)適合的模型的殘差序列應(yīng)是白噪聲序列?？蓱?yīng)用Box-LjungQ統(tǒng)計(jì)量對殘差進(jìn)行檢驗(yàn)，檢驗(yàn)水準(zhǔn)α=0.05。

⑤模型的預(yù)測：利用所建模型進(jìn)行預(yù)測，并與實(shí)際數(shù)據(jù)進(jìn)行比較，評(píng)價(jià)模型的優(yōu)劣。

2 結(jié)果

2.1序列的平穩(wěn)化圖1是2004年1月至2011年12月我國丙肝月發(fā)病數(shù)(Yt)的時(shí)間序列圖。從圖1可以看出，該序列按年呈現(xiàn)出明顯的非平穩(wěn)性和季節(jié)性，隨著時(shí)間t的增大，Yt的波動(dòng)越來越大，即呈現(xiàn)遞增型的異方差。丙肝發(fā)病數(shù)顯現(xiàn)總體上升趨勢，發(fā)病數(shù)在每年的3、4月份達(dá)到高峰，隨后慢慢下降，在次年1月份左右達(dá)到低谷。

為消除異方差，對原始數(shù)據(jù)進(jìn)行對數(shù)變換，變換后丙肝月發(fā)病數(shù)據(jù)(lnYt)的ACF和PACF圖見圖2?？梢钥闯?，lnYt是非平穩(wěn)的，且存在周期為12個(gè)月的季節(jié)波動(dòng)。因此對lnYt進(jìn)行一階非季節(jié)差分和一階季節(jié)差分，得到新的序列ΔΔ12lnYt(差分算子Δ=1-L,Δ12=1-L12)。從ΔΔ12lnYt的ACF和PACF圖(圖3)可以看出，該序列近似為一個(gè)平穩(wěn)過程。

2.2模型的識(shí)別與定階由于原始序列對數(shù)變換后，經(jīng)過一階非季節(jié)差分和一階季節(jié)差分達(dá)到平穩(wěn)，因此s=12,d=1,D=1，模型可初步確定為ARIMA(p,1,q)(P,1,Q)12，其中p、q、P和Q待定。觀察序列ΔΔ12lnYt的PACF圖，其偏相關(guān)函數(shù)和自相關(guān)函數(shù)均在滯后1階以后降為0，因此p=1，q=1。

2.3參數(shù)估計(jì)及檢驗(yàn)通過計(jì)算，9個(gè)模型中備選模型ARIMA(1,1,1)(1,1,0)12和ARIMA(1,1,1)(2,1,0)12的參數(shù)具有統(tǒng)計(jì)學(xué)意義，結(jié)果見表1。其中AR(1)表示階數(shù)為1的非季節(jié)自回歸過程，SAR(12)表示周期為12、階數(shù)為1的季節(jié)自回歸過程，SAR(24)表示周期為12、階數(shù)為2的季節(jié)自回歸過程，MA(1)表示階數(shù)為1的非季節(jié)移動(dòng)平均過程。另外，模型所有根(包括實(shí)根和復(fù)根)的倒數(shù)均小于1，均符合建模要求。

圖1 丙肝月發(fā)病數(shù)(Yt)的時(shí)間序列圖

圖2 lnYt的PACF(上)和ACF(下)圖

圖3 ΔΔ12lnYt的PACF(上)和ACF(下)圖

從表1可以看出，模型ARIMA(1,1,1)(2,1,0)12的R2和調(diào)整R2均較ARIMA(1,1,1)(1,1,0)12大，而統(tǒng)計(jì)量AIC比ARIMA(1,1,1)(1,1,0)12小，提示模型ARIMA(1,1,1)(2,1,0)12的擬合效果較好。

根據(jù)表1估計(jì)的參數(shù)值得到模型的表達(dá)式為[3-4]：(1+0.222L)(1+0.820L12+0.694L24)(1-L)(1-L12)lnYt=(1+0.648L)εt,其中L為后移算子，LkYt=Yt-k，εt為白噪聲序列。

2.4模型的診斷在Eviews 6.0中對模型ARIMA(1,1,1)(2,1,0)12的殘差進(jìn)行Q檢驗(yàn)， 結(jié)果P均>0.05，說明模型ARIMA(1,1,1)(2,1,0)12的殘差序列是白噪聲序列，提示所選模型恰當(dāng)，適用于預(yù)測。

2.5模型的預(yù)測用模型ARIMA(1,1,1)(2,1,0)12對原序列進(jìn)行擬合，結(jié)果見圖4。由圖4可以看出，預(yù)測值與實(shí)際值基本吻合，動(dòng)態(tài)趨勢基本一致。用模型ARIMA(1,1,1)(2,1,0)12對我國2012年1月至12月丙肝逐月發(fā)病人數(shù)進(jìn)行預(yù)測，并用實(shí)際數(shù)據(jù)進(jìn)行預(yù)測精度的驗(yàn)證，結(jié)果見表2。由表2可以看出，該模型的預(yù)測結(jié)果相對誤差最大為0.176，最小為0.009，平均為0.078，預(yù)測效果較好。

2.6發(fā)病規(guī)律的周期性動(dòng)態(tài)變化趨勢與模型預(yù)測效果的關(guān)系丙肝和甲肝2004年1月至2011年12月發(fā)病數(shù)歸一化后的方差計(jì)算結(jié)果見表3。丙肝的平均方差為0.030，甲肝的平均方差為0.047，丙肝的方差更小，所以可以認(rèn)為丙肝原始數(shù)據(jù)周期性動(dòng)態(tài)變化比甲肝更趨一致。用模型ARIMA(1,1,1)(2,1,0)12對我國2012年1月至12月甲肝逐月發(fā)病人數(shù)的預(yù)測結(jié)果相對誤差最大為0.266，最小為0.003，平均為0.138，與表2結(jié)果比較，可以看出丙肝的預(yù)測效果明顯優(yōu)于甲肝，平均相對誤差減小了近一半。

圖4 實(shí)際值與預(yù)測值序列比較圖

表2 用模型ARIMA(1,1,1)(2,1,0)12對我國2012年1月至12月丙肝逐月發(fā)病人數(shù)的預(yù)測結(jié)果

表3 丙肝和甲肝2004年1月至2011年12月發(fā)病數(shù)歸一化后各月的方差

3 討論

近年來，應(yīng)用時(shí)間序列模型對傳染病進(jìn)行監(jiān)控、預(yù)測得到了較廣泛的應(yīng)用。在傳染病的預(yù)測研究方面，目前用的較多的模型有微分方程模型、灰色預(yù)測模型、Markov模型、通徑分析模型等[7-11]。微分方程模型是比較簡單和基本的封閉體系模型，它適合描述傳染病的自然發(fā)展過程，但在現(xiàn)實(shí)生活中人們會(huì)通過藥物、隔離治療等方式人為地控制傳染病的發(fā)展，這會(huì)在一定程度上影響預(yù)測效果?；疑A(yù)測模型是在灰色理論的基礎(chǔ)上對生成數(shù)列建模，只需很少數(shù)據(jù)量即可完成預(yù)測，但當(dāng)數(shù)據(jù)波動(dòng)比較大時(shí)，預(yù)測精確度就大大下降。Markov模型是利用概率建立一種隨機(jī)時(shí)序模型來進(jìn)行預(yù)測，預(yù)測準(zhǔn)確度依賴于概率的準(zhǔn)確性以及時(shí)間序列的獨(dú)立性。通徑分析模型是回歸分析的補(bǔ)充和發(fā)展，需要結(jié)合專業(yè)知識(shí)建造通徑圖，當(dāng)變量較多時(shí)，何種通徑圖較優(yōu)則很難確定。以上模型都無法對周期性數(shù)據(jù)進(jìn)行建模和預(yù)測。

ARIMA模型比較靈活，具有不受數(shù)據(jù)類型束縛和適用性強(qiáng)的特點(diǎn)，既適用于非周期性序列，也適用于周期性序列，且周期可以為年份、季度、月份，故適用范圍廣泛。ARIMA季節(jié)模型整合了趨勢因素、周期因素和隨機(jī)誤差等因素的原始時(shí)間序列變量，通過差分?jǐn)?shù)據(jù)轉(zhuǎn)換等方法將非平穩(wěn)序列轉(zhuǎn)變?yōu)榱憔档钠椒€(wěn)隨機(jī)序列，通過反復(fù)識(shí)別和模型診斷、比較，選擇理想的模型進(jìn)行數(shù)據(jù)擬合和預(yù)測。該方法結(jié)合了自回歸和移動(dòng)平均方法的長處，可不受數(shù)據(jù)類型的束縛，適用性較強(qiáng)，是一種短期預(yù)測效果較好的模型，非常適用于與季節(jié)周期性相關(guān)疾病的研究。傳染病的發(fā)病往往與季節(jié)有關(guān)，具有較明顯的周期性。作者采用ARIMA季節(jié)模型對我國丙肝的疫情動(dòng)態(tài)和發(fā)展趨勢進(jìn)行預(yù)測，得到了較好的預(yù)測效果，這為我國丙肝發(fā)病的預(yù)測提供了一種可行的方法。

如何更好地應(yīng)用ARIMA季節(jié)模型對傳染病發(fā)病狀況進(jìn)行預(yù)測，作者對發(fā)病規(guī)律的周期性動(dòng)態(tài)變化趨勢對模型預(yù)測效果的影響進(jìn)行了分析。作者對我國2004年1月至2011年12月丙肝和甲肝的發(fā)病數(shù)據(jù)按年歸一化處理后計(jì)算方差，結(jié)果丙肝的方差更小，說明丙肝原始數(shù)據(jù)周期性動(dòng)態(tài)變化比甲肝更趨一致。用ARIMA季節(jié)模型對2012年1月至12月逐月發(fā)病人數(shù)進(jìn)行預(yù)測，丙肝預(yù)測結(jié)果的平均相對誤差較甲肝減小了近一半，說明丙肝的預(yù)測效果明顯優(yōu)于甲肝。該研究結(jié)果證實(shí)ARIMA季節(jié)模型的預(yù)測效果與樣本數(shù)據(jù)的周期性動(dòng)態(tài)變化趨勢有關(guān)，周期性變化趨勢越一致，預(yù)測結(jié)果也越準(zhǔn)確。這為今后更好地開展工作、提高預(yù)測的精度提供了科學(xué)的依據(jù)。

[1]中華預(yù)防醫(yī)學(xué)會(huì)醫(yī)院感染控制分會(huì).中國丙型病毒性肝炎醫(yī)院感染防控指南(2012年)[J].中國醫(yī)學(xué)前沿雜志：電子版,2012，4(11)：57

[2]Bowerman BL, O'Connell RT.Forecasting and time series:an applied approach[M].3rd ed.北京：機(jī)械工業(yè)出版社，2003：437

[3]方積乾，陸盈，張晉昕，等.現(xiàn)代醫(yī)學(xué)統(tǒng)計(jì)學(xué)[M].北京：人民衛(wèi)生出版社，2002.

[4]Box GEP,Jenkins GM.Time series analysis： forecasting and control[M].Holden-Day,1976.

[5]攸頻，張曉峒.Eviews 6實(shí)用教程[M].北京：中國財(cái)政經(jīng)濟(jì)出版社，2008.

[6]易丹輝.數(shù)據(jù)分析與EViews應(yīng)用[M].北京:中國人民大學(xué)出版社，2008.

[7]張彥琦，唐責(zé)立，王文昌，等.ARIMA模型及其在肺結(jié)核預(yù)測中的應(yīng)用[J].現(xiàn)代預(yù)防醫(yī)學(xué)，2008，35(9)：1608

[8]胡建利，祖榮強(qiáng)，彭志行，等.江蘇省戊型肝炎發(fā)病趨勢的時(shí)間序列模型應(yīng)用[J].南京醫(yī)科大學(xué)學(xué)報(bào)：自然科學(xué)版，2011，31(12)：1874

[9]胡建利,梁祁,吳瑩,等.季節(jié)時(shí)間序列模型在菌痢發(fā)病預(yù)測中的應(yīng)用[J].中國衛(wèi)生統(tǒng)計(jì)，2012，29(1)：34

[10]Muhammad HL,Nur Haizum Suhartono AR,Mohd TL,et al.Seasonal ARIMA for forecasting air pollution index:a case study[J].Am J Applied Sci，2012，9(4)：570

[11]王丙剛,曲波,郭海強(qiáng),等.傳染病預(yù)測的數(shù)學(xué)模型研究[J].中國衛(wèi)生統(tǒng)計(jì),2007,24(5):536

(2013-07-23收稿 責(zé)任編輯王 曼)

Application of seasonal ARIMA model in forecasting incidence of hepatitis C in China

YULinfeng1),WUJing2),ZHOUSuolan1),DINGYong2)

1)DepartmentofBiomedicalEngineering,NanjingMedicalUniversity,Nanjing210029 2)DepartmentofMathematicsandComputerSciences,NanjingMedicalUniversity,Nanjing210029

seasonal ARIMA model; hepatitis C; incidence; prediction

Aim: To forecast the incidence of hepatitis C in China using seasonal ARIMA model.Methods: Seasonal ARIMA model was established based on the monthly reported cases data of hepatitis C in China from 2004 to 2011,and used to forecast the data of 2012.Actual data of 2012 were used to assess prediction effect. The model establishment and forecasting for hepatitis A were carried out using the same method. The variance of hepatitis A and hepatitis C incidence from 2004 to 2011 normalized according to the years was calculated. The predicted effect of hepatitis A and hepatitis C was compared.Results: The model of ARIMA(1,1,1)(2,1,0)12was established successfully.The expression of the model was (1+0.222L)(1+0.820L12+0.694L24)(1-L)(1-L12)lnYt=(1+0.648L)εt,the parameters were as follows: AR(1)=-0.222(t=-2.392,P=0.020)，SAR(12)=-0.820(t=-8.009,P<0.001),SAR(24)=-0.694(t=-6.124,P<0.001),MA(1)= -0.648(t=-5.889,P<0.001),residual error sequence was white noise sequence (P>0.05), theR2of fitting was 0.824 and the averge error of prediction was 0.078. The averge variances of hepatitis C and hepatitis A normalized incidence were 0.030 and 0.047, suggesting that periodic dynamic change of hepatitis C data was more consistent. The averge relative error of prediction of hepatitis A was 0.138,higher than that of hepatitis C.Conclusion: ARIMA(1,1,1)(2,1,0)12season model can be used to predict incidence of hepatitis C in China. Periodic dynamic change trend of sample data is more consistent, the ARIMA seasonal model predicted result is more accurate.

10.13705/j.issn.1671-6825.2014.03.014

*江蘇省大學(xué)生實(shí)踐創(chuàng)新訓(xùn)練計(jì)劃項(xiàng)目 2012JSSPITP1033；南京醫(yī)科大學(xué)基礎(chǔ)醫(yī)學(xué)院優(yōu)勢學(xué)科教師培養(yǎng)基金項(xiàng)目 JX10131801099

R512.6