2004-2018年我國狂犬病疫情時間序列分析

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

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

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

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

1 材料和方法

1.1資料來源 2004-2018年全國狂犬病月發病統計數據來自中國疾病預防控制中心“疾病監測信息報告管理系統”。

1.2 分析方法

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

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

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

1.2.4模型的驗證和評價 以我國2018年1-12月狂犬病月發病數據為驗證樣本，以平均絕對誤差和平均相對誤差為評價標準，將最優模型所得的預測結果和實際結果進行比較，評價最優模型的預測精準度。

1.2.5模型的應用 利用最優模型對我國2019年狂犬病疫情進行預測。

2 結 果

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

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

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

圖2 2005-2017年我國狂犬病月發病數據經對數轉換、一階普通差分和一階季節性差分后序列圖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年我國狂犬病月發病數呈現出明顯的季節性變化趨勢，故采用ARIMA乘積季節模型對我國狂犬病疫情進行擬合建模。對原始時間序列進行一階普通差分和一階季節性差分后獲得平穩時間序列，故可確定d=1，D=1，ARIMA乘積季節模型表現為：ARIMA(p,1,q)×(P,1,Q)12。一般情況下，p、q、P、Q不超過二階，故對所有符合條件的模型進行篩選。此外，根據對數轉換和一階差分處理后的自相關和偏相關分析圖，以正態化的BIC、平均絕對標準化誤差、均方根誤差和標準化的R2為參考依據，初步篩選出5個擬合度較高的模型，其結果見表1。

圖3a: 差分處理前自相關系數分析圖；圖3b: 差分處理前偏自相關系數分析圖；圖3c: 差分處理后自相關系數分析圖；圖3d: 差分處理后偏自相關系數圖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 差分處理前后自相關系數和偏自相關系數分析圖Fig.3 Autocorrelation coefficient (ACF) and the partial autocorrelation coefficient (PACF) of monthly rabies incidence time series before (after) differential processing

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

Tab.1 Fitting parameters of different ARIMA models

ARIMA模型模型擬合指標平穩化的R2RMSEMAPEMAE標準化的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(決定系數)，RMSE(均方誤差平方根)，MAPE(平均絕對誤差百分比)，MAE(平均絕對誤差)

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

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

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

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

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

2.4 模型的驗證和評價

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

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

表3 2018年我國狂犬病發病預測值和實際值

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

月份實際值預測值絕對誤差相對誤差%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模型的應用 利用最優模型ARIMA(0,1,1)×(0,1,1)12對我國2019年狂犬病月發病數據進行預測，結果顯示，2019年我國狂犬病總發病數預計達208例。

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

3 討 論

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

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

利益沖突：無