ARIMA乘積季節(jié)模型在青州市布魯氏菌病發(fā)病預(yù)測(cè)中的應(yīng)用

劉 杰 武欽發(fā) 李偉國 肖宇飛 毛 倩 石福艷 王素珍

1 山東省青州市地方病防治研究所，262500 山東 青州； 2 山東省青州市疾病預(yù)防控制中心，262500 山東 青州；3 濰坊醫(yī)學(xué)院衛(wèi)生統(tǒng)計(jì)學(xué)教研室，261053 山東 濰坊

布魯氏菌病是由布魯氏菌屬細(xì)菌入侵機(jī)體引起人和牛、羊、豬等動(dòng)物共患的傳染病[1]。人布魯氏菌病被《中華人民共和國傳染病防治法》列入乙類法定報(bào)告?zhèn)魅静」芾?。青州市首次?965年發(fā)現(xiàn)布魯氏菌病病人，其后布魯氏菌病在青州市的流行基本處于穩(wěn)定狀態(tài)，但是2011年以來,由于養(yǎng)殖業(yè)的不斷發(fā)展、牲畜交易市場(chǎng)開放、相應(yīng)的防治措施落實(shí)不到位，布魯氏菌病疫情有不斷上升趨勢(shì)。為了更好地掌握青州市布魯氏菌病流行趨勢(shì)，本文根據(jù)青州市2011—2017年布魯氏菌病月發(fā)病數(shù)建立自回歸求和移動(dòng)平均(ARIMA)乘積季節(jié)模型，并應(yīng)用該模型預(yù)測(cè)青州市2018年各月份布魯氏菌病發(fā)病數(shù)，用實(shí)際月發(fā)病數(shù)進(jìn)行回代擬合，檢驗(yàn)?zāi)Ｐ偷念A(yù)測(cè)效果,以探討ARIMA 模型預(yù)測(cè)布魯氏菌病月發(fā)病數(shù)的可行性。

1 資料與方法

1.1 資料來源

青州市2011—2017年布魯氏菌病疫情資料來源于中國疾病預(yù)防控制系統(tǒng)，人口學(xué)資料來源于青州市統(tǒng)計(jì)年鑒和2011—2017年青州市國民經(jīng)濟(jì)和社會(huì)發(fā)展公報(bào)。

1.2 ARIMA模型簡介

1.2.1定義

ARIMA模型全稱為自回歸移動(dòng)平均模型，記作ARIMA(p,d,q)，由西方學(xué)者George Edward Pelham Box和 Gwilym Meirion Jenkins于20世紀(jì)70年代提出，是一種著名的時(shí)間序列統(tǒng)計(jì)模型[2]。其含義是指將非平穩(wěn)時(shí)間序列轉(zhuǎn)化為平穩(wěn)時(shí)間序列，然后將因變量僅對(duì)它的滯后值以及隨機(jī)誤差項(xiàng)的現(xiàn)值和滯后值進(jìn)行回歸所建立的模型。

1.2.2建模步驟

第1步，時(shí)間序列平穩(wěn)化。首先繪制時(shí)間序列的散點(diǎn)圖或折線圖，粗略判斷其平穩(wěn)性。對(duì)于非平穩(wěn)的時(shí)間序列，如果存在異方差，則應(yīng)對(duì)序列進(jìn)行對(duì)數(shù)轉(zhuǎn)換；如果時(shí)間序列存在一定變化趨勢(shì)，則應(yīng)對(duì)序列進(jìn)行差分處理，差分包括季節(jié)性差分和非季節(jié)性差分，一般取d與D值為1或2，使其轉(zhuǎn)為平穩(wěn)時(shí)間序列[2-3]。轉(zhuǎn)換后的時(shí)間序列是否平穩(wěn)可通過單位根檢驗(yàn)(augmented dickey-fuller test, ADF)來判定。

第2步，模型識(shí)別與定階。模型的識(shí)別就是根據(jù)時(shí)間序列的識(shí)別規(guī)則判斷模型適用于哪種過程，建立擬合效果最優(yōu)的模型。用到的工具為自相關(guān)函數(shù)ACF和偏自相關(guān)函數(shù)PACF以及它們各自的相關(guān)圖。通過觀察分析自相關(guān)函數(shù)和偏自相關(guān)函數(shù)的拖尾性和結(jié)尾性判斷時(shí)間序列適用哪種模型[3-4]。對(duì)于季節(jié)性ARIMA模型，先識(shí)別非季節(jié)成分，確定p、q值，再識(shí)別季節(jié)成分，P、Q階數(shù)較難確定，一般來說超過2 階的情況很少見。可將階數(shù)從低階到高階進(jìn)行不同的組合，選擇出幾個(gè)粗模型[3,5]。

第3步，參數(shù)估計(jì)與診斷。在得到幾個(gè)備選模型后，需進(jìn)行參數(shù)估計(jì)，檢驗(yàn)是否有統(tǒng)計(jì)學(xué)意義。參數(shù)估計(jì)的方法有極大似然估計(jì)、最小二乘法估計(jì)、矩估計(jì)等統(tǒng)計(jì)方法[2]。經(jīng)過參數(shù)估計(jì)得到幾個(gè)備選模型后，用Ljung-Box 統(tǒng)計(jì)量對(duì)模型進(jìn)行適應(yīng)性檢驗(yàn)，檢查該模型時(shí)間序列所蘊(yùn)含的信息提取是否充分，若P>0.05則該模型殘差是白噪聲。同時(shí)根據(jù)AIC 和BIC 準(zhǔn)則判定模型的擬合優(yōu)度。經(jīng)過反復(fù)試驗(yàn)比較，將模型系數(shù)有統(tǒng)計(jì)學(xué)意義、殘差序列是白噪聲、AIC 和BIC值較小的模型定為預(yù)測(cè)模型[5-6]。

第4步，預(yù)測(cè)及應(yīng)用。利用選擇出的最優(yōu)模型進(jìn)行預(yù)測(cè)，將預(yù)測(cè)值與實(shí)際值進(jìn)行比較，評(píng)價(jià)模型的預(yù)測(cè)效果。通過計(jì)算平均絕對(duì)百分比誤差(MAPE)、均方根誤差(RMSE)等指標(biāo)來評(píng)估模型預(yù)測(cè)的準(zhǔn)確性和精度[7-8]。

1.3 統(tǒng)計(jì)學(xué)方法

本文采用Excel 2013 對(duì)2011—2017年青州市布魯氏菌病月發(fā)病數(shù)據(jù)進(jìn)行整理，采用SPSS 25.0軟件構(gòu)建青州市布魯氏菌病的ARIMA 模型。α=0.05。

2 結(jié)果

2.1 序列平穩(wěn)化

本研究以青州市2011—2017年布魯氏菌病月發(fā)病數(shù)為基礎(chǔ)數(shù)據(jù)，用2018年1—12月份青州市布魯氏菌病實(shí)際月發(fā)病數(shù)對(duì)構(gòu)建的ARIMA模型進(jìn)行驗(yàn)證。首先繪制時(shí)間序列圖(圖1)，圖1顯示青州市布魯氏菌病發(fā)病數(shù)總體呈現(xiàn)上升趨勢(shì)，且具備季節(jié)性和周期性特點(diǎn)，發(fā)病集中在4—7月份，5月份為發(fā)病高峰，提示該時(shí)間序列是不平穩(wěn)的時(shí)間序列，需要對(duì)該時(shí)間序列平穩(wěn)化處理。由于存在周期性，對(duì)該序列同時(shí)進(jìn)行一般性差分和季節(jié)性差分，周期為12，即d和D的取值為1。經(jīng)差分后的時(shí)間序列圖見圖2，此時(shí)序列已基本平穩(wěn)。

2.2 模型識(shí)別和參數(shù)估計(jì)

對(duì)獲得的平穩(wěn)時(shí)間序列繪制自相關(guān)函數(shù)圖ACF和偏自相關(guān)函數(shù)圖PACF(見圖3)。由ACF圖和PACF圖可見，從一階開始后均視為拖尾，因此建立ARIMA(p,1,q)(P,1,Q)12模型，其中p和q取值為0或1；根據(jù)經(jīng)驗(yàn)，P和Q超過2階的情況很少，采取從低階到高階依次組合的方法逐一進(jìn)行嘗試，應(yīng)用Box-Ljung 方法對(duì)殘差序列進(jìn)行白噪聲檢驗(yàn)剔除非白噪聲模型 (P<0.05) ，初步篩選出4 種平穩(wěn)的白噪聲模型ARIMA 模型。這4種模型的標(biāo)準(zhǔn)化BIC值和Box-Ljung 檢驗(yàn)結(jié)果見表1。經(jīng)過對(duì)比，這4種模型中ARIMA(0,1,1)(1,1,0)12模型BIC值相對(duì)較小，固定R平方相對(duì)較大并且該模型參數(shù)均具有統(tǒng)計(jì)學(xué)意義(P<0.01),該模型參數(shù)估計(jì)見表2。確定模型ARIMA(0,1,1)(1,1,0)12為最優(yōu)模型。繪制ARIMA(0,1,1)(1,1,0)12模型的殘差序列 ACF 和PACF 圖(圖4)對(duì)模型進(jìn)行診斷，本模型殘差序列 ACF 和 PACF 圖顯示殘差基本都落在 95%區(qū)間內(nèi)(P>0.05)。表明該模型的信息提取充分，殘差為白噪聲。初步認(rèn)為該模型包含原始時(shí)間序列的所有特征，并且時(shí)間序列之間不存在相關(guān)性，可以用來預(yù)測(cè)青州市2018年各個(gè)月份人布魯氏菌病月發(fā)病報(bào)告數(shù)。

圖1 青州市2011—2017年布魯氏菌病月發(fā)病人數(shù)序列圖

圖2 經(jīng)過差分后青州市2011—2017年布魯氏菌病月發(fā)病人數(shù)序列

圖3 經(jīng)過差分轉(zhuǎn)換后的青州市布魯氏菌病月發(fā)病數(shù)自相關(guān)圖和偏自相關(guān)圖

表1 4種備選ARIMA模型的BIC值和Ljung-Box Q結(jié)果

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

圖4 殘差序列的自相關(guān)(ACF)和偏自相關(guān)(PACF)圖

2.3 模型預(yù)測(cè)及應(yīng)用

采用青州市2011年1月至2017年12月布魯氏菌病月發(fā)病人數(shù)回代入ARIMA(0,1,1)(1,1,0)12模型進(jìn)行擬合，繪制擬合曲線，見圖5。結(jié)果顯示擬合值與觀測(cè)值總體趨勢(shì)一致，擬合值與觀測(cè)值重合程度高，模型的MAPE值為59.23%，觀測(cè)值和擬合值都在預(yù)測(cè)值的95%置信區(qū)間內(nèi)，絕對(duì)誤差和相對(duì)誤差均較小，說明該模型的預(yù)測(cè)精度較高，可以對(duì)未來進(jìn)行很好地跟蹤和預(yù)測(cè)。用該模型對(duì)青州市2018年各月布魯氏菌病的新發(fā)病人數(shù)進(jìn)行預(yù)測(cè)，并與實(shí)際月發(fā)病數(shù)進(jìn)行比較，結(jié)果見表3。表3顯示，實(shí)際值與預(yù)測(cè)值的平均誤差為1.1，預(yù)測(cè)趨勢(shì)與實(shí)際趨勢(shì)基本一致，模型的預(yù)測(cè)效果較好。

圖5 2011—2017年青州市布魯氏菌病實(shí)際值與模型預(yù)測(cè)值比較

表3 2018年青州市布魯氏菌病實(shí)際月發(fā)病數(shù)與模型預(yù)測(cè)月發(fā)病數(shù)比較

3 討論

布魯氏菌病在青州市的流行可追溯到1965年[9-10]，對(duì)青州市人群的健康影響較大。ARIMA模型是一種著名的時(shí)間序列統(tǒng)計(jì)模型,不需要分析影響疾病的各種因素，將各種影響因素的綜合效應(yīng)統(tǒng)一蘊(yùn)涵于時(shí)間變量之中，研究時(shí)間變量的歷史數(shù)據(jù)，挖掘出其中的規(guī)律以準(zhǔn)確預(yù)測(cè)未來的發(fā)展變化趨勢(shì)的數(shù)理統(tǒng)計(jì)模型[11]。

本文采用2011—2017年青州市布魯氏菌病發(fā)病資料經(jīng)過序列平穩(wěn)化、模型識(shí)別和參數(shù)估計(jì)等一系列過程構(gòu)建了ARIMA(0,1,1)(1,1,0)12預(yù)測(cè)模型，得到2018年布魯氏菌病各月發(fā)病數(shù)預(yù)測(cè)值，并與實(shí)際值比較。預(yù)測(cè)結(jié)果顯示布魯氏菌病預(yù)測(cè)值和實(shí)際值的動(dòng)態(tài)趨勢(shì)基本一致，實(shí)際值均落入預(yù)測(cè)值的95%置信區(qū)間內(nèi),表明ARIMA(0,1,1)(1,1,0)12模型擬合效果較好，精度較高，體現(xiàn)出模型良好的實(shí)用性和應(yīng)用價(jià)值[12]。因此ARIMA模型可以為傳染病發(fā)病率進(jìn)行早期預(yù)警，為傳染病防控工作提供科學(xué)依據(jù)，具體可根據(jù)傳染病既往的發(fā)展變化規(guī)律判斷暴發(fā)或者流行的可能性，如果實(shí)際值在預(yù)測(cè)值95%置信區(qū)間范圍內(nèi)波動(dòng)，表明疫情基本正常;如果超出預(yù)測(cè)值95%置信區(qū)間, 提示傳染病暴發(fā)或流行的可能性較大[13]。

盡管ARIMA模型兼具回歸分析和移動(dòng)平均的優(yōu)點(diǎn)，適用范圍較因果回歸分析法等分析方法廣，研究過程相對(duì)簡化、經(jīng)濟(jì)、適用，但是ARIMA模型是假定時(shí)間序列未來發(fā)展模式是由過去的慣性趨勢(shì)發(fā)展來的，因此只適用于短期預(yù)測(cè)。而且ARIMA模型對(duì)數(shù)據(jù)的要求較高，要求時(shí)間序列在平穩(wěn)的前提下要有30個(gè)以上數(shù)據(jù)。另外當(dāng)實(shí)際情況復(fù)雜，特別是在采取了干預(yù)措施(如加強(qiáng)管理)或外部環(huán)境等因素發(fā)生較大改變時(shí)，模型的建立相對(duì)比較困難，應(yīng)根據(jù)實(shí)際情況綜合考慮預(yù)測(cè)結(jié)果[14]?？梢栽谒Ｐ椭胁粩嘣黾訕颖玖繉?duì)模型進(jìn)行修正甚至重新擬合，進(jìn)一步提高模型的預(yù)測(cè)精度。

綜上所述，本研究建立的ARIMA(0,1,1)(1,1,0)12模型能反映青州市布魯氏菌病時(shí)間分布特征和發(fā)展趨勢(shì)，預(yù)測(cè)效果較好，可以為疫情監(jiān)測(cè)和疾病防控工作提供參考依據(jù)。