謝 淵，劉淑清，董國英，朱武洋

狂犬病是一種由狂犬病病毒引起，以犬、狼、貓等食肉動物為主要傳播媒介，以恐水、怕風(fēng)、進行性癱瘓為主要臨床特征的急性和致命性的人獸共患病[1-2]。狂犬病的癥狀一旦發(fā)展，其死亡率幾乎為100%，目前該病仍高居甲乙類傳染病死亡率首位[3]。狂犬病主要通過動物的咬傷或刮傷進行傳播，狂犬病病毒感染中樞神經(jīng)系統(tǒng)最終導(dǎo)致腦部疾病和死亡[4-5]。所有種類的哺乳動物對狂犬病病毒易感，但狗仍是狂犬病的主要載體且大部分人間狂犬病都由其引起[6-7]。

狂犬病在全球范圍內(nèi)流行，每年導(dǎo)致數(shù)萬人死亡，而95%以上發(fā)生在亞洲和非洲[8]，其中印度和中國是報道病例數(shù)最多的國家[9]。公元前556年我國首次出現(xiàn)狂犬病病例報道，目前該病仍是非常嚴重的公共衛(wèi)生問題。為進一步了解我國狂犬病疫情的分布特征和流行趨勢，并對狂犬病疫情進行短期內(nèi)預(yù)測，現(xiàn)利用2004-2018年我國狂犬病發(fā)病數(shù)據(jù)建立季節(jié)性時間序列并對其進行分析。

自回歸移動平均模型(Autoregressive Integrated Moving Average Model, ARIMA) 即Box-Jenkins模型，一般表現(xiàn)形式為：ARIMA(p,d,q)×(P,D,Q)s。其中p(P)、d(D)、q(Q)分別為非季節(jié)性(季節(jié)性)的自回歸平均階數(shù)、移動平均階數(shù)、差分次數(shù)，s則為模型的季節(jié)周期[10]。該模型在時間序列分析中應(yīng)用廣泛，對傳染病的短期預(yù)測具有良好的擬合效果[11-12]。故本研究基于ARIMA模型，對2004-2018年我國狂犬病疫情進行時間序列分析，并對其進行短期預(yù)測以期對我國狂犬病的防控提供參考。

1 材料和方法

1.1資料來源 2004-2018年全國狂犬病月發(fā)病統(tǒng)計數(shù)據(jù)來自中國疾病預(yù)防控制中心“疾病監(jiān)測信息報告管理系統(tǒng)”。

1.2 分析方法

1.2.1序列的建立和平穩(wěn)化 將2004-2017年我國狂犬病月發(fā)病數(shù)據(jù)時間單位定義為年份、季度和月份，而后可得到相應(yīng)的時間序列曲線圖，通過時間序列圖觀察序列的平穩(wěn)性，利用SPSS 19.0對不穩(wěn)定的時間序列數(shù)據(jù)進行數(shù)據(jù)轉(zhuǎn)化和差分處理達到序列平穩(wěn)化的目的[13]，并達到以下要求：均數(shù)和方差不隨時間變化；自相關(guān)系數(shù)僅與時間間隔相關(guān)。

1.2.2模型的識別和定階 狂犬病疫情呈現(xiàn)季節(jié)性變化特征，故可用ARIMA模型進行擬合。ARIMA模型一般表現(xiàn)形式為：ARIMA(p,d,q)×(P,D,Q)s。其中p(P)、d(D)、q(Q)分別為非季節(jié)性(季節(jié)性)的自回歸平均階數(shù)、移動平均階數(shù)、差分次數(shù)，s則為模型的季節(jié)周期。定階過程即根據(jù)平穩(wěn)時間序列的自相關(guān)圖(ACF)和偏自相關(guān)圖(PACF)初步確定模型參數(shù)的過程。通過觀察ACF和PACF的截尾或拖尾的情況對模型進行擬合，比較所得到的擬合結(jié)果并對其做出相應(yīng)的調(diào)整，初步建立一個或多個可擬合的ARIMA模型。

1.2.3模型的檢驗和優(yōu)化 根據(jù)平穩(wěn)的R2、正態(tài)化的BIC和平均絕對標(biāo)準(zhǔn)化誤差(MASE)等指標(biāo)對模型的擬合度進行檢測評價。同時對模型進行Ljung-Box檢驗，判斷模型殘差序列是否為白噪聲序列，篩選出通過檢驗的模型，確定正態(tài)化的BIC值最小的為最優(yōu)模型。

1.2.4模型的驗證和評價 以我國2018年1-12月狂犬病月發(fā)病數(shù)據(jù)為驗證樣本，以平均絕對誤差和平均相對誤差為評價標(biāo)準(zhǔn)，將最優(yōu)模型所得的預(yù)測結(jié)果和實際結(jié)果進行比較，評價最優(yōu)模型的預(yù)測精準(zhǔn)度。

1.2.5模型的應(yīng)用 利用最優(yōu)模型對我國2019年狂犬病疫情進行預(yù)測。

2 結(jié) 果

2.1序列的建立和平穩(wěn)化 以2004-2017年我國狂犬病月發(fā)病數(shù)據(jù)建立時間序列圖(圖1)，橫軸表示2004年1月-2017年12月時間軸，縱軸表示期間每年月發(fā)病數(shù)。所建立的時間序列圖顯示2004-2007年我國狂犬病發(fā)病呈現(xiàn)上升趨勢，至2007年達到高峰，當(dāng)年全國報告發(fā)病人數(shù)達3 300人，而2008-2017年我國狂犬病發(fā)病數(shù)逐漸下降。圖中還顯示每年的8-10月發(fā)病數(shù)可達至高峰，說明8-10月為狂犬病的高發(fā)月份?？傮w而言，2004年1月-2017年12月間，我國狂犬病發(fā)病數(shù)波動較大，呈現(xiàn)出以年為周期的變化趨勢，整體疫情存在明顯的季節(jié)性變化，由此說明該時間序列為非平穩(wěn)時間序列。

觀察原始序列圖發(fā)現(xiàn)序列呈現(xiàn)周期性變化規(guī)律，且周期為12個月。因此為了獲得平穩(wěn)的時間序列，需對原始序列進行對數(shù)轉(zhuǎn)換和差分處理。通過嘗試后發(fā)現(xiàn)，經(jīng)自然對數(shù)轉(zhuǎn)換、一階普通差分和一階季節(jié)性差分后可消除原始序列的變化趨勢從而使之達到平穩(wěn)的狀態(tài)，得到的序列圖如圖2所示。圖2表明差分后時間序列圖的長期趨勢和季節(jié)性變化趨勢消失，且數(shù)值在零上下波動。此外，觀察差分前后的自相關(guān)系數(shù)和偏自相關(guān)系數(shù)分析圖(圖3)發(fā)現(xiàn)，僅當(dāng)k=1和k=12時，自相關(guān)系數(shù)突破了置信區(qū)間，自相關(guān)系數(shù)和偏自相關(guān)系數(shù)在k=1后逐漸呈現(xiàn)衰減的趨勢且逐漸落入可信區(qū)間范圍內(nèi)，故認為此時的時間序列已趨于平穩(wěn)。

圖1 2004-2017我國狂犬病月發(fā)病數(shù)據(jù)時間序列圖Fig.1 Time series of monthly incidence data of rabies in China from 2004 to 2017

圖2 2005-2017年我國狂犬病月發(fā)病數(shù)據(jù)經(jīng)對數(shù)轉(zhuǎn)換、一階普通差分和一階季節(jié)性差分后序列圖Fig.2 Time series of monthly incidence data of rabies after logarithmic transformation, first-order ordinary difference and first-order seasonal difference, 2005-2017

2.2模型的識別和定階 由于2004-2017年我國狂犬病月發(fā)病數(shù)呈現(xiàn)出明顯的季節(jié)性變化趨勢，故采用ARIMA乘積季節(jié)模型對我國狂犬病疫情進行擬合建模。對原始時間序列進行一階普通差分和一階季節(jié)性差分后獲得平穩(wěn)時間序列，故可確定d=1，D=1，ARIMA乘積季節(jié)模型表現(xiàn)為：ARIMA(p,1,q)×(P,1,Q)12。一般情況下，p、q、P、Q不超過二階，故對所有符合條件的模型進行篩選。此外，根據(jù)對數(shù)轉(zhuǎn)換和一階差分處理后的自相關(guān)和偏相關(guān)分析圖，以正態(tài)化的BIC、平均絕對標(biāo)準(zhǔn)化誤差、均方根誤差和標(biāo)準(zhǔn)化的R2為參考依據(jù)，初步篩選出5個擬合度較高的模型，其結(jié)果見表1。

圖3a: 差分處理前自相關(guān)系數(shù)分析圖；圖3b: 差分處理前偏自相關(guān)系數(shù)分析圖；圖3c: 差分處理后自相關(guān)系數(shù)分析圖；圖3d: 差分處理后偏自相關(guān)系數(shù)圖Fig.3a: Autocorrelation coefficient (ACF) of monthly rabies incidence time series; Fig.3b: Partial autocorrelation coefficient (PACF) of monthly rabies incidence time series; Fig.3c: Autocorrelation coefficient (ACF) of monthly rabies incidence time series after differential processing; Fig.3d: Partial autocorrelation coefficient (PACF) of monthly rabies incidence time series after differential processing 圖3 差分處理前后自相關(guān)系數(shù)和偏自相關(guān)系數(shù)分析圖Fig.3 Autocorrelation coefficient (ACF) and the partial autocorrelation coefficient (PACF) of monthly rabies incidence time series before (after) differential processing

表1 不同自回歸階數(shù)和移動平均階數(shù)ARIMA模型的擬合參數(shù)

Tab.1 Fitting parameters of different ARIMA models

ARIMA模型模型擬合指標(biāo)平穩(wěn)化的R2RMSEMAPEMAE標(biāo)準(zhǔn)化的BICARIMA(0,1,1)×(0,1,1)120.58819.77810.02414.2076.164ARIMA(0,1,1)×(2,1,1)120.58820.12810.69014.3096.265ARIMA(0,1,2)×(0,1,1)120.58819.83510.66914.2806.202ARIMA(0,1,1)×(1,1,1)120.58620.12710.67614.3076.232ARIMA(0,1,1)×(1,1,0)120.46322.42812.14815.9676.416

注:R2(決定系數(shù))，RMSE(均方誤差平方根)，MAPE(平均絕對誤差百分比)，MAE(平均絕對誤差)

2.3模型的檢驗和優(yōu)化 正態(tài)化的BIC越小，標(biāo)準(zhǔn)化的R2越大，模型擬合效果越好[14]，因此可以確定擬合效果最好的模型為ARIMA(0,1,1)×(0,1,1)12。利用Box-Ljung方法對所獲得的最優(yōu)模型的殘差序列進行白噪聲檢驗，檢驗結(jié)果顯示最優(yōu)模型Ljung-Box Q=14.413，自由度為16，P>0.05，無統(tǒng)計學(xué)意義，表明該模型的殘差序列為白噪聲序列。ARIMA(0,1,1)×(0,1,1)12模型的參數(shù)估計結(jié)果(表2)顯示，MA滯后和MA季節(jié)性滯后的估計值均有統(tǒng)計學(xué)意義(P<0.05)。此外，根據(jù)最優(yōu)模型的殘差A(yù)CF和PACF分析圖(圖4)可知，殘差A(yù)CF和PACF大部分都處于95%置信區(qū)間內(nèi)，表明該殘差的分布是隨機的，不存在相關(guān)性，滿足獨立性檢驗。綜上所述，獲得的ARIMA(0,1,1)×(0,1,1)12模型是有效的且擬合效果較好。

表2 模型ARIMA(0,1,1)×(0,1,1)12的參數(shù)估計結(jié)果

Tab.2 Parameter estimation result of ARIMA(0,1,1)×(0,1,1)12

參數(shù)類別遲滯模型系數(shù)t值P值常數(shù)0.2900.935-0.3100.757MA滯后10.6810.06410.5950.000MA,季節(jié)性滯后10.9940.0910.4750.003季節(jié)性差分1

圖4 ARIMA(0,1,1)×(0,1,1)12模型殘差序列的自相關(guān)系數(shù)和偏自相關(guān)系數(shù)圖Fig.4 Autocorrelation coefficient and partial autocorrelation coefficient of ARIMA(0,1,1)×(0,1,1)12

2.4 模型的驗證和評價

2.4.1回代擬合 將獲得的最優(yōu)模型ARIMA(0,1,1)×(0,1,1)12對2004-2017年的月發(fā)病數(shù)據(jù)進行回代擬合(圖5)，結(jié)果顯示擬合值和觀測值變化趨勢基本一致且觀測值一直處于預(yù)測可信區(qū)間內(nèi)。

2.4.2模型預(yù)測 利用獲得的最優(yōu)模型ARIMA(0,1,1)×(0,1,1)12對2018年我國狂犬病月發(fā)病數(shù)據(jù)進行預(yù)測，將預(yù)測值與真實值進行比較，結(jié)果如表3所示。結(jié)果表明預(yù)測結(jié)果的平均相對誤差為14.91%，按預(yù)測12個月總發(fā)病數(shù)據(jù)來看，相對誤差為0.82%。

表3 2018年我國狂犬病發(fā)病預(yù)測值和實際值

Tab.3 Predictive and actual values of rabies in China in 2018

月份實際值預(yù)測值絕對誤差相對誤差%1月2833517.862月2924-517.243月1927842.114月2832414.295月3033310.006月343400.007月1821316.678月4137-49.769月383800.0010月4236-614.2811月3225-721.8712月2723-414.81總計366363-30.82

2.5模型的應(yīng)用 利用最優(yōu)模型ARIMA(0,1,1)×(0,1,1)12對我國2019年狂犬病月發(fā)病數(shù)據(jù)進行預(yù)測，結(jié)果顯示，2019年我國狂犬病總發(fā)病數(shù)預(yù)計達208例。

圖5 模型ARIMA(0,1,1)×(0,1,1)12回代擬合結(jié)果Fig.5 Fitting result of ARIMA(0,1,1)×(0,1,1)12

3 討 論

狂犬病是一種在全球范圍內(nèi)廣泛流行的急性傳染病和人獸共患病，全球已有150余個國家深受其害，每年造成大量的經(jīng)濟損失。截至目前，亞洲和非洲仍然是狂犬病高發(fā)地區(qū)，狂犬病仍在人間、家養(yǎng)動物和野生動物之間循環(huán)[15]。我國是繼非洲后狂犬病疫情最為嚴重的國家，自2007年以來，我國狂犬病疫情呈現(xiàn)緩解的趨勢，但形勢依舊嚴峻。因此，狂犬病的疫情監(jiān)測和防控仍是我國傳染病防制工作的重點之一。

時間序列分析是一種定量預(yù)測的方法，它可將疾病發(fā)生的多種影響因素如地理環(huán)境和季節(jié)環(huán)境等綜合考慮于時間變量中，分析發(fā)病數(shù)據(jù)隨時間變化的規(guī)律進而對疾病的發(fā)生進行短期內(nèi)的預(yù)測[16]。目前，ARIMA模型是時間序列分析最基礎(chǔ)、最常見的方法，該模型不僅可對狂犬病疫情進行預(yù)測，也廣泛應(yīng)用于其他傳染病的預(yù)測分析[17-19]。本研究采用了ARIMA乘積季節(jié)模型對我國2004-2018年狂犬病月發(fā)病情況進行了時間序列分析，目的在于分析近年來我國狂犬病的流行特征和流行趨勢，并對狂犬病的短期流行做出預(yù)測。分析結(jié)果顯示模型預(yù)測的我國狂犬病發(fā)病趨勢和實際值吻合度較高，相對誤差較小，表明ARIMA模型可用于我國狂犬病疫情的預(yù)測。但由于ARIMA模型是利用歷史數(shù)據(jù)進行統(tǒng)計分析，且疾病的發(fā)生還受其他諸多因素的影響，故ARIMA模型只適用于疾病的短期預(yù)測[20]。因此，為了對我國狂犬病的疫情進行更加精準(zhǔn)的預(yù)測，需不斷補充新的數(shù)據(jù)，結(jié)合狂犬病發(fā)生的其他影響因素，調(diào)整模型的參數(shù)以適應(yīng)狂犬病的實際發(fā)生情況。

利益沖突：無