2004-2018年我國(guó)狂犬病疫情時(shí)間序列分析

謝 淵，劉淑清，董國(guó)英，朱武洋

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

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

自回歸移動(dòng)平均模型(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ù)、移動(dòng)平均階數(shù)、差分次數(shù)，s則為模型的季節(jié)周期[10]。該模型在時(shí)間序列分析中應(yīng)用廣泛，對(duì)傳染病的短期預(yù)測(cè)具有良好的擬合效果[11-12]。故本研究基于ARIMA模型，對(duì)2004-2018年我國(guó)狂犬病疫情進(jìn)行時(shí)間序列分析，并對(duì)其進(jìn)行短期預(yù)測(cè)以期對(duì)我國(guó)狂犬病的防控提供參考。

1 材料和方法

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

1.2 分析方法

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

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

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

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

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

2 結(jié) 果

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

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

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

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

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

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

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

參數(shù)類(lèi)別遲滯模型系數(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 模型的驗(yàn)證和評(píng)價(jià)

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

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

表3 2018年我國(guó)狂犬病發(fā)病預(yù)測(cè)值和實(shí)際值

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

月份實(shí)際值預(yù)測(cè)值絕對(duì)誤差相對(duì)誤差%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總計(jì)366363-30.82

2.5模型的應(yīng)用 利用最優(yōu)模型ARIMA(0,1,1)×(0,1,1)12對(duì)我國(guó)2019年狂犬病月發(fā)病數(shù)據(jù)進(jìn)行預(yù)測(cè)，結(jié)果顯示，2019年我國(guó)狂犬病總發(fā)病數(shù)預(yù)計(jì)達(dá)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余個(gè)國(guó)家深受其害，每年造成大量的經(jīng)濟(jì)損失。截至目前，亞洲和非洲仍然是狂犬病高發(fā)地區(qū)，狂犬病仍在人間、家養(yǎng)動(dòng)物和野生動(dòng)物之間循環(huán)[15]。我國(guó)是繼非洲后狂犬病疫情最為嚴(yán)重的國(guó)家，自2007年以來(lái)，我國(guó)狂犬病疫情呈現(xiàn)緩解的趨勢(shì)，但形勢(shì)依舊嚴(yán)峻。因此，狂犬病的疫情監(jiān)測(cè)和防控仍是我國(guó)傳染病防制工作的重點(diǎn)之一。

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

利益沖突：無(wú)