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

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

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

資料與方法

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

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

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

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

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

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

結(jié) 果

1.一般情況

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

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

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

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

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

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

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

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

表2 備選模型及殘差白噪聲檢驗(yàn)表

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

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

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

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

討 論

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

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

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