ARIMA季節(jié)性模型在預(yù)測蘭州市丙肝發(fā)病人數(shù)中的應(yīng)用*

馬 晶 王 梅 譙小偉 曹文珮 李娟生 任曉衛(wèi) 王宇紅 劉小寧△

【提 要】 目的 應(yīng)用自回歸移動平均(autoregressive integrated moving average，ARIMA)乘積季節(jié)模型建立適合預(yù)測蘭州市丙肝月報告發(fā)病人數(shù)的最優(yōu)模型，為丙肝防控工作提供參考依據(jù)。方法 收集2010年1月至2019年12月蘭州市丙肝月報告發(fā)病人數(shù)，基于2010年1月至2018年12月丙肝月報告發(fā)病人數(shù)建立ARIMA乘積季節(jié)模型，比較2019年1-12月預(yù)測值和實際值以進行預(yù)測模型性能檢驗，并對2020年進行預(yù)測。結(jié)果 蘭州市丙肝月報告發(fā)病人數(shù)總體呈現(xiàn)下降趨勢，并于每年2月份達到最低值，于3月份迅速上升至峰值，有季節(jié)性和周期性。根據(jù)2010年1月至2018年12月丙肝月報告發(fā)病數(shù)據(jù)擬合出最佳預(yù)測模型ARIMA(2，1，0)(0，1，1)12，該模型擬合2019年1-12月的預(yù)測值與實際值的相對誤差范圍是0.00% ～ 57.28%，平均相對誤差為16.96%，2020年預(yù)測結(jié)果顯示蘭州市丙肝報告發(fā)病人數(shù)略有下降。結(jié)論 基于本研究數(shù)據(jù)，ARIMA(2，1，0)(0，1，1)12模型能較好地擬合蘭州市丙肝的月報告發(fā)病人數(shù)，可用于短期的預(yù)測。

丙型肝炎(簡稱丙肝)是由丙型肝炎病毒(hepatitis C virus，HCV)引起的一種以肝臟損害為主的傳染性疾病，該病毒可造成急性或慢性肝炎感染，其嚴(yán)重程度從持續(xù)幾周的輕微病癥到終身嚴(yán)重疾病不等[1]。2017年世界衛(wèi)生組織(WHO)估計，全球HCV感染率為1.0%，約有7100萬慢性HCV感染者，HCV已成為當(dāng)今世界嚴(yán)重的公共衛(wèi)生問題之一[2]。近年來，中國HCV報告發(fā)病數(shù)呈逐年上升趨勢，全國報告發(fā)病率由2010年的11.41/10萬上升至2019年的16.02/10萬，HCV的防控形勢依然嚴(yán)峻[3]。在目前無安全有效疫苗的情況下，應(yīng)充分利用傳染病監(jiān)測信息資源，開展有效的HCV發(fā)病狀況預(yù)測，這對HCV防控工作的開展具有重要的指導(dǎo)意義。本研究擬采用自回歸移動平均(auto-regressive integrated moving average，ARIMA)乘積季節(jié)模型建立蘭州市HCV月報告發(fā)病數(shù)模型，在此基礎(chǔ)上進行擬合和短期預(yù)測，為有關(guān)部門制定HCV防治策略提供科學(xué)的參考依據(jù)。

資料與方法

1.資料來源 2010年1月至2019年12月蘭州市HCV監(jiān)測數(shù)據(jù)來源于中國疾病預(yù)防控制信息系統(tǒng)平臺下的傳染病信息直報系統(tǒng)。

2.診斷標(biāo)準(zhǔn) HCV報告病例診斷依據(jù)《丙型肝炎診斷標(biāo)準(zhǔn)(WS213-2018)》[4]。

3.質(zhì)量控制 納入標(biāo)準(zhǔn)：現(xiàn)住址為蘭州市、已通過審核的臨床診斷病例和確診病例(病例均為單純感染HCV)，剔除信息缺失病例。

4.研究方法 ARIMA模型是時間序列分析中最受歡迎和方便的模型之一，并且廣泛用于傳染病的預(yù)測中[5]。本研究選用ARIMA(p，d，q)(P，D，Q)s模型建立預(yù)測蘭州市HCV月報告發(fā)病數(shù)的數(shù)學(xué)模型，其中p、d、q分別為自回歸階數(shù)、非季節(jié)差分階數(shù)、移動平均階數(shù)，P、D、Q分別代表季節(jié)自回歸階數(shù)、季節(jié)差分階數(shù)、季節(jié)性移動平均階數(shù)，s為周期。

ARIMA模型建模的基本步驟[6-8]：(1)時間序列的平穩(wěn)化。模型要求原始序列平穩(wěn)，即均值和方差不隨時間變化，因此通過對原始時間序列作時序圖、自相關(guān)圖(autocorrelation function，ACF)、偏自相關(guān)圖(partial autocorrelation function，PACF)以及單位根檢驗(augmented dickey-fuller test，ADF)判斷原始時間序列是否穩(wěn)定。若原序列不平穩(wěn)，則需要對非平穩(wěn)序列通過對數(shù)轉(zhuǎn)換、趨勢差分、周期性差分等預(yù)處理轉(zhuǎn)變?yōu)闊o長期趨勢性和周期性的平穩(wěn)序列。(2)模型識別。對預(yù)處理后達到平穩(wěn)性要求的序列，繪制ACF圖和PACF圖，根據(jù)其截尾或拖尾的情況識別模型中p、q、P、Q的數(shù)值，擬合一個或多個ARIMA模型。(3)參數(shù)估計和模型診斷。運用最大似然估計法(maximum likelihood estimation，MLE)對模型參數(shù)進行估計，對備選模型進行Ljung-Box殘差白噪聲檢驗(L-B檢驗)，以判斷模型中的信息提取是否充分，再通過比對備選模型的池赤信息準(zhǔn)則(Akaike information criterion，AIC)值和施瓦茨貝葉斯準(zhǔn)則(Schwarz′s Bayesian criterion，SBC)值，在所有通過檢驗的模型中，AIC或SBC函數(shù)達到最小的模型為相對最優(yōu)模型。(4)模型預(yù)測及評價。利用選出的最優(yōu)模型進行預(yù)測，將預(yù)測值和實際值進行比較，使用相對誤差評價模型的準(zhǔn)確性。

5.統(tǒng)計分析 采用excel 2019建立蘭州市HCV月報告發(fā)病人數(shù)數(shù)據(jù)庫，運用R 4.0.1軟件進行ARIMA模型構(gòu)建、分析及預(yù)測。其中2010年1月至2018年12月HCV月報告發(fā)病數(shù)作為訓(xùn)練集擬合ARIMA乘積季節(jié)模型，2019年1-12月的數(shù)據(jù)作為測試集檢驗ARIMA模型的外推能力。以P<0.05為差異有統(tǒng)計學(xué)意義。

結(jié) 果

1.一般情況

2010-2019年蘭州市HCV累計報告病例數(shù)22695例，年平均報告人數(shù)2269例，年均報告發(fā)病率為70.20/10萬，見表1。對2010年1月-2018年12月HCV月報告發(fā)病數(shù)的時間序列進行趨勢性、季節(jié)性和隨機誤差分解，結(jié)果顯示，HCV月報告發(fā)病數(shù)在2010-2012年呈平緩上升趨勢，2013年出現(xiàn)下降趨勢，說明此時間序列為非平穩(wěn)時間序列。在季節(jié)因素的作用下，HCV月報告發(fā)病例數(shù)在每年的2月份出現(xiàn)低谷，隨后3月份出現(xiàn)高峰，5、8、12月份報告人數(shù)有增加的趨勢，隨機誤差維持在一定范圍內(nèi)。見圖1。

表1 2010-2019年蘭州市HCV報告發(fā)病情況

圖1 原始數(shù)據(jù)的時間序列圖

2.ARIMA模型構(gòu)建

(1)數(shù)據(jù)預(yù)處理 由原始序列圖可以看出此時間序列為非平穩(wěn)序列，且具有一定的季節(jié)性。使用“diff”語句對原始序列進行1階12步季節(jié)差分，消除時間序列趨勢和季節(jié)影響后，使用tseries包中的“ adf.test ”語句進行ADF檢驗(P=0.01)，由此可知差分后數(shù)據(jù)滿足平穩(wěn)性要求。使用“Box.test”語句對差分處理后的數(shù)據(jù)進行白噪聲檢驗，序列延遲6階、12階的檢驗結(jié)果分別為χ2=50.92(P<0.001)、χ2=94.19(P<0.001)，其為非白噪聲序列，說明差分處理后的序列值之間還蘊含著較強的相關(guān)關(guān)系，必須進一步分析。因此，可考慮建立ARIMA(p，1，q)(P，1，Q)s模型進行分析和預(yù)測。

(2)模型識別 ARIMA模型中參數(shù)p和q的識別主要依據(jù)ACF圖和PACF圖的截尾或拖尾情況，對差分后的數(shù)據(jù)分別做ACF圖(圖2b)與PACF圖(圖2c)，ACF圖顯示自相關(guān)系數(shù)0階截尾(q=0)，且在12階、24階均有明顯波動，提示差分后序列有季節(jié)效應(yīng)。PACF圖顯示偏自相關(guān)系數(shù)在2階后落入置信區(qū)間，即p=2。根據(jù)差分次數(shù)和數(shù)據(jù)的周期情況初步判定模型為ARIMA(2，1，0)(P，1，Q)12，采用從低階向高階的順序進行嘗試，根據(jù)模型的擬合優(yōu)度指標(biāo)選擇最佳模型。

圖2 1階12步差分變換后分析圖

(3)參數(shù)估計和模型診斷 運用 R 4.0.1軟件中的“auto.arima()”函數(shù)自動識別出的模型為ARIMA(0，1，1)(0，0，2)12，基于以上分析共有10個備選模型，通過計算，10個備選模型中有6個模型參數(shù)具有統(tǒng)計學(xué)意義，見表2。比較發(fā)現(xiàn)， ARIMA(2，1，0)(0，1，1)12模型的AIC和SBC值最小，故選取該模型作為本次研究相對最優(yōu)模型。經(jīng)檢驗，ARIMA(2，1，0)(0，1，1)12模型各參數(shù)檢驗統(tǒng)計量P值均小于0.01，L-B檢驗結(jié)果顯示序列延遲6、12、18及24階差異均無統(tǒng)計學(xué)意義(P>0.05)，說明殘差序列為白噪聲序列，模型擬合較好，見表3。

表2 備選模型及殘差白噪聲檢驗表

表3 ARIMA(2，1，0)×(0，1，1)12模型參數(shù)和殘差白噪聲檢驗結(jié)果

(4)模型預(yù)測效果評估 為驗證ARIMA(2，1，0)(0，1，1)12模型的外推能力，用“forecast”函數(shù)預(yù)測2019年1-12月蘭州市HCV的月報告發(fā)病例數(shù)，預(yù)測結(jié)果見表4。結(jié)果顯示，預(yù)測HCV月發(fā)病例數(shù)相對誤差范圍為0.00%～57.28%，平均相對誤差為16.96%。2019年1-12月中除9月份外其他月份HCV報告發(fā)病的實際值均落在預(yù)測值的95%置信區(qū)間(95%CI)范圍內(nèi)，提示該模型預(yù)測性能較好。對2020年HCV月發(fā)病數(shù)進行預(yù)測(圖3)，結(jié)果顯示2020年蘭州市HCV月報告發(fā)病數(shù)略有下降。

表4 2019年1-12月HCV的月發(fā)病例數(shù)的實際值與預(yù)測值比較

圖3 蘭州市HCV月報告發(fā)病人數(shù)預(yù)測

討 論

本研究結(jié)果顯示，蘭州市HCV流行主要有以下特征：第一，蘭州市HCV報告發(fā)病人數(shù)總體呈下降趨勢。從原始時序分解圖中觀察到蘭州市HCV報告發(fā)病人數(shù)在2010-2012年呈現(xiàn)上升趨勢，2013-2014年HCV報告發(fā)病人數(shù)較2012年明顯下降，提示可能與開展HCV數(shù)據(jù)質(zhì)量核查、規(guī)范HCV病例報告等相關(guān)，也可能與同期甘肅省衛(wèi)生廳出臺規(guī)范乙肝診斷與報告的文件間接相關(guān)[9-10]。2014-2019年HCV報告發(fā)病人數(shù)下降趨勢逐漸平緩，由此可見，隨著HCV防治宣傳教育的普及，人群關(guān)于HCV知識的知曉率不斷提升，防范意識也逐漸增強。第二，存在季節(jié)性升高。通過本研究原始序列圖可觀察到HCV發(fā)病周期為一年，季節(jié)性特征明顯，可觀察到除2012年外蘭州市每年HCV的月報告發(fā)病人數(shù)在2月份均處于較低水平，其后迅速上升，在3月份會達到峰值，形成“春節(jié)效應(yīng)”[11-12]，提示可能與HCV病例因我國春節(jié)影響審核且上報延遲以及春節(jié)過后人口流動量增加有關(guān)。本研究結(jié)果還顯示，每年5月份、8月份和12月份報告人數(shù)也會出現(xiàn)增加的趨勢，同樣提示可能有周期性或季節(jié)性的存在。但也有研究指出HCV的感染不具有季節(jié)性[13-15]，提示可能與HCV主要為血源及性傳播疾病有關(guān)，但季節(jié)性不明顯不代表其發(fā)病無高峰和低谷期，因此仍要做好相關(guān)的防控工作。

結(jié)果顯示，2019年1-12月HCV報告發(fā)病的實際值基本落在預(yù)測值的95%CI內(nèi)，從預(yù)測精度來看，該模型預(yù)測的最小相對誤差為0，最大相對誤差為57.28%，平均相對誤差為16.96%。雖然該模型能較好地擬合HCV發(fā)病人數(shù)，但部分月份如9月份的擬合值和實際值誤差較大，提示可能與該月HCV報告發(fā)病人數(shù)變化過大有關(guān)。HCV的感染通常還涉及吸毒、艾滋病以及性傳播方式等，因此HCV的報告發(fā)病人數(shù)不僅受檢測技術(shù)、疾病報告系統(tǒng)完善性等客觀條件的影響，還受到患者行為方式、文化程度以及宗教信仰等因素的影響[16]，而ARIMA模型的優(yōu)勢是將多種因素的綜合效應(yīng)蘊含于時間變量之中，在不考慮這些因素的情況下對疾病的發(fā)生情況仍能進行科學(xué)的預(yù)測，并通過反復(fù)識別和修正達到最滿意的預(yù)測效果[17]。此外，目前HCV的治療仍缺乏特效方法，也沒有有效的疫苗進行預(yù)防，且HCV是一種慢性化程度高、危害極為嚴(yán)重的傳染病，感染后極易轉(zhuǎn)化為慢性肝炎、肝硬化、肝癌等經(jīng)濟負擔(dān)較嚴(yán)重的疾病。因此，通過利用現(xiàn)有的監(jiān)測數(shù)據(jù)，觀察各個時間點上的變化規(guī)律，并對該類疾病在發(fā)病高峰前進行積極主動的預(yù)測，根據(jù)預(yù)測發(fā)出預(yù)警，這對疾病防控具有重要的意義。

ARIMA模型在疾病發(fā)病預(yù)測中也存在一定的偏倚。首先，時間序列是按照一定歷史時期的規(guī)律進行預(yù)測，其前提是這些規(guī)律在一定時期內(nèi)保持相對穩(wěn)定，但是當(dāng)未來某個時點或時段的實際值變化過大時，模型的預(yù)測值可能就會偏離實際值，故預(yù)測時間不宜太長。其次，在模型識別時有主觀成分存在，由于個人經(jīng)驗的不同，不同的研究者對同一組數(shù)據(jù)建模時可能得不到相同的最優(yōu)模型，這也會導(dǎo)致偏倚的發(fā)生，從而影響模型預(yù)測的精確度和準(zhǔn)確度。由于HCV感染是一個復(fù)雜的過程，因此，在今后的研究中應(yīng)該多考慮其他因素對預(yù)測結(jié)果的影響，建立預(yù)測效果更好、更適合于HCV預(yù)測的模型。