全國人間狂犬病疫情的時間序列分析

，，，，，　，　，　

狂犬病是狂犬病病毒感染引起的人獸共患病。除南極洲外全球各大陸均有病例報告，每年約有6萬例狂犬病死亡，其中超過95%的死亡在亞洲和非洲[1]。WHO估計約99%的人狂犬病感染來源為犬。此病目前無特異療法，病死率幾乎100%。為了解我國狂犬病疫情時間分布特征及流行趨勢，現用全國疫情資料進行分析研究。

ARIMA模型即Box-Jenkins模型。一般形式：ARIMA(p,d,q)×(P,D,Q)s，其中d 和D分別為差分和季節性差分次數；p和P分別是自回歸(AR)模型和季節性AR模型的參數個數；q和Q分別是移動平均 (MA)模型和季節性MA模型的參數個數；S是季節性期數。對有季節性變化的時間序列應建立ARIMA(p,d,q)× (P,D,Q)s 模型。

1　資料與方法

1.1資料來源2004年1月至2016年11月全國狂犬病月發病數來自國家衛計委(http:// www.nhfpc.gov.cn/zhuz/index.shtml)，2004-2015年全國人口數據來自國家統計局(http://www.stats.gov.cn/)。

1.2分析方法① 序列建立及特征分析：據月發病數及人口數計算月發病率，并建立2004年1月至2016年12月發病率時間序列。用LOESS光滑對序列做季節性分解，同時計算各月的季節指數，據季節指數確定發病高峰。② 序列的平穩化：用自相關圖、偏自相關圖及單位根檢驗對序列進行平穩性檢驗；對非平穩序列用差分法平穩化。平穩序列用Liung-Box檢驗是否為白噪聲；對白噪聲序列選擇合適的分析方法；對非白噪聲序列建立模型。③ 模型的識別及建立：據平穩時間序列的自相關圖(ACF)和偏自相關圖(PACF)，初步估計模型階數。按AIC值最小原則選擇最優模型。④ 模型擬合及評價：根據建立的模型擬合模型參數、繪制殘差的QQ圖，檢驗是否滿足正態分布。用Liung-Box檢驗殘差是否滿足獨立性要求。⑤ 預測：用平均絕對標準化誤差(mean absolute error,MASE)衡量模型的擬合優度，MASE<1時模型擬合較好。用建立的模型預測2016年1-11月全國狂犬病月發病率，計算預測值和實際值間的相對誤差以判斷模型的預測精度。

2　結　果

2.1疫情概況2004年1月至2016年11月的月發病率時間序列(圖1A)。用LOESS光滑對序列進行季節性分解結果顯示，疫情存在明顯的季節性，其中8～10月為高峰(32.4%)，季節指數分別為1.277、1.285、1.224(圖1B)；自2007年后發病率呈逐年下降趨勢(圖1C)。

圖1-A：月發病率時間分布；圖1-B：季節效應；圖1-C：年度趨勢分布；圖1-D：殘差隨機波動Fig.1-A：temporal distribution graph；Fig.1-B：seasonality effect graph；Fig.1-C：annual trend；Fig.1-D：residual distribution圖1　全國狂犬病月發病率序列季節性分解Fig.1　The seasonality decomposition diagram of monthly rabies incidence time series of China

2.2模型識別及建立用2004-2015年月發病率數據建立模型。其自相關圖(ACF)顯示為非平穩序列，需進行差分(圖2)。對序列進行1階差分后，顯示為平穩序列(Dickey-Fuller=-9.352，P<0.01)；進行Box-Liung檢驗顯示為非白噪聲序列(C2=4.906,P=0.0268)。由于有季節性高峰，故用ARIMA(p,d,q)(P,D,Q)12模型,因此d=1，D=0。據1階差分序列的ACF圖和PACF圖的圖形類別，逐步選擇合適的p、q、P及D值，以AIC值最小者為最優模型，結果顯示最優模型為ARIMA(2,1,1)(2,0,0)12，AIC=-57.07。

圖2　全國狂犬病2004-2015年月發病率時間序列自相關圖ACF (左)和偏自相關圖PACF(右)Fig.2　The autocorrelation function (ACF,left) and the partial autocorrelation function (PACF,right) of the monthly rabies incidences in China (2004-2015)

2.3模型擬合及評價對建立的ARIMA(2,1,1)(2,0,0)12進行擬合，所有系數均有統計學意義(表1)。繪制模型殘差的QQ圖顯示滿足正態分布(圖3)。模型殘差ACF圖及PACF圖顯示無明顯的相關性(圖4)；Box-Liung檢驗結果也顯示模型殘差無明顯的相關性(C2=0.339,P=0.560)、是白噪聲序列。

表1ARIMA(2,1,1)(2,0,0)12模型參數擬合值
Tab.1The fitting values of parameters in ARIMA(2,1,1)(2,0,0)12

　　指標indexAR1aAR2bMA1cSAR1dSAR2e系數coefficient0．6650．202-0．9850．4520．372標準誤standarderror0．0960．0980．0310．0800．086t值tvalue6．927??2．061?3．177??5．65??4．326??

注釋：a和b分別為AR模型的第1和第2個參數，c為MA模型的第1個參數，d和e分別是季節性AR模型的第1和第2個參數；*P<0.05；**P<0.01。

圖3　ARIMA(2,1,1)(2,0,0)12模型的殘差QQ圖Fig.3　Residual QQ plot of ARIMA(2,1,1)(2,0,0)12

圖4　模型殘差的自相關圖ACF(左)和偏自相關圖PACF(右)Fig.4　The autocorrelation function (ACF,left) and the partial autocorrelation function (PACF,right) of the model residual

2.4預測 R軟件分析結果顯示模型ARIMA(2,1,1)(2,0,0)12的MASE=0.755,此值<1表明擬合效果較好。用此模型預測2016年1-11月發病率，結果相對誤差-17.25%～62.48%，平均相對誤差為15.61%(表2)。

表2全國狂犬病2016年1-11月份月發病率預測結果
Tab.2The predicted incidences of rabies in China,from Jan.to Nov.,2016

月份Month實際值actualvalue預測值(95%CI)predictedvalue(95%CI)絕對誤差absoluteerror相對誤差relativeerror10．4000．511(0．146～0．875)0．1110．27620．3350．355(0～0．795)0．01990．05930．2330．378(0～0．882)0．1450．62540．2110．309(0～0．858)0．0980．46550．4730．391(0～0．975)-0．082-0．17360．3930．501(0～1．112)0．1080．27470．4290．498(0～1．132)0．0690．16080．4580．570(0～1．222)0．1120．24490．5530．543(0～1．211)-0．010-0．018100．5600．537(0～1．217)-0．023-0．041110．3780．521(0～1．212)0．1430．377

3　討　論

全球除印度外中國年報告狂犬病死亡最多[1]。50年代以來，我國人間狂犬病出現3次高峰[2-3]。第1個高峰50年代中期年死亡最高達1 900多例。第2次高峰為80年代初期，1981年死亡7 037例是最高年份，80年代年均死亡5 537例。第3次高峰在21世紀初期，疫情重現快速增長趨勢。主要宿主動物是犬，調查顯示93.74%病例的感染來源是犬[4]。廣西、貴州、廣東、湖南和四川是近年來高發省份。此病潛伏期數天至≥1年，但多為20～90 d[2]。可在每年高峰季節到來前的3～4個月，對<15歲兒童、農民及流動人口等高危人群及高危地區，加強健康教育，提高暴露后的規范處置率，以減少發病。

時間序列分析據該觀察值歷史數據來預測其發展趨勢，探索其時間分布特征及影響因素。ARIMA模型應用于經濟、工程、氣象、水利等眾多領域。醫學衛生領域也用于預測疾病發病率、人群壽命、醫療衛生費用和食物中毒等，還用于研究氣候及環境因素對人群健康的影響[5-11]。

本研究用2004年1月至2016年11月全國狂犬病月發病率數據進行時間序列分析。顯示2007年以來發病率呈下降趨勢，疫情有明顯季節性(8～10月高峰)，與既往報道的時間分布特征類似[2,4]。

建立ARIMA模型至少要有50 個時間點或7～8 個周期的數據。本研究用了144個數據。通過數據差分、模型識別，建立最優模型ARIMA(2,1,1)(2,0,0)12。已知該模型的MASE=0.755，擬合效果較好。對2016年1-11月發病率預測顯示效果較好，平均相對誤差15.61%。由于ARIMA模型預測時是假定事件影響因素無明顯變化，因此ARIMA模型常用于時間序列數據的短期預測。

參考文獻：

[1] Jackson AC.Human rabies: a 2016 update[J].Curr Infect Dis Rep,2016,18(11):38.DOI: 10.1007/s11908-016-0540-y

[2] Zhou H,Vong S,Liu K,et al.Human rabies in China,1960-2014: a descriptive epidemiological study[J].PLoS Negl Trop Dis,2016,10(8):e4874.DOI: 10.1371/journal.pntd.0004874

[3] Mao WC,Lan RW,Meng NX,et al.Rabies epidemical and prevention and control research progress[J].J Med Pest Control,2015(3):283-286.(in Chinese)

毛偉成,藍榮偉,蒙南新,等.狂犬病流行與防治研究進展[J].醫學動物防制,2015(3):283-286.

[4] Song M,Tang Q,Rayner S,et al.Human rabies surveillance and control in China,2005-2012[J].BMC Infect Dis,2014,14:212.DOI:10.1186/1471-2334-14-212

[5] Wang KW,Deng C,Li JP,et al.Hybrid methodology for tuberculosis incidence time-series forecasting based on ARIMA and a NAR neural network[J].Epidemiol Infect,2017:1-12.DOI:10.1017/S0950268816003216

[6] 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-151.DOI:10.1017/S0950268815001144

[7] Anwar MY,Lewnard JA,Parikh S,et al.Time series analysis of malaria in Afghanistan: using ARIMA models to predict future trends in incidence[J].Malar J,2016,15(1):566.DOI:10.1186/s12936-016-1602-1

[8] Zhang X,Pang Y,Cui M,et al.Forecasting mortality of road traffic injuries in China using seasonal autoregressive integrated moving average model[J].Ann Epidemiol,2015,25(2):101-106.DOI:10.1016/j.annepidem.2014.10.015

[9] Chadsuthi S,Iamsirithaworn S,Triampo W,et al.Modeling seasonal influenza transmission and its association with climate factors in Thailand using time-series and ARIMAX analyses[J].Comput Math Methods Med,2015,2015:436495.DOI: 10.1155/2015/436495

[10] Zhang Y,Wang SN,Liu Y,et al.Application of ARIMA model on predicting monthly hospital admissions and hospitalization expenses for respiratory diseases[J].Chin J Health Statistic,2015，32(2):197-200.(in Chinese)

張越，王勝難，劉媛，等.應用ARIMA模型對呼吸系統疾病月住院量及住院費用的預測[J].中國衛生統計,2015，32(2):197-200.

[11] Wang L,Yao WQ.Forecasting the epidemic trend of hand-foot-and-mouth disease with time sequence model in Liaoning.[J].Chin J Health Statistic,2016，33(5):847-849.(in Chinese)

王伶，姚文清.利用時間序列模型分析預測遼寧手足口病疫情趨勢[J].中國衛生統計,2016，33(5):847-849.