遼寧省疾病預(yù)防控制中心(110005) 常 琳 孫 靜 方 興

【提 要】 目的 建立疑似預(yù)防接種異常反應(yīng)(adverse event following immunization，AEFI)報(bào)告病例數(shù)的季節(jié)性自回歸積分滑動(dòng)平均模型(autoregressive integrated moving average model，ARIMA)，對(duì)AEFI報(bào)告數(shù)進(jìn)行預(yù)測(cè)。方法 使用R語言做模型的識(shí)別、模型參數(shù)的估計(jì)、檢驗(yàn)，建立季節(jié)性ARIMA模型，對(duì)遼寧省AEFI報(bào)告數(shù)進(jìn)行模型擬合，用2019年1-12月的預(yù)測(cè)值與實(shí)際值作比較，檢驗(yàn)?zāi)Ｐ偷念A(yù)測(cè)能力。結(jié)果 經(jīng)過多次檢驗(yàn)，確定ARIMA(0，1，1)(1，1，1)12模型預(yù)測(cè)能力最佳，其殘差序列為白噪聲。用2019年1-12月數(shù)據(jù)檢驗(yàn)?zāi)Ｐ?，由MAPE的絕對(duì)值可以看出，除3月外其他月份預(yù)測(cè)值與實(shí)際值相差均較小，說明模型的擬合優(yōu)度相對(duì)較好，預(yù)測(cè)結(jié)果可靠。結(jié)論 季節(jié)性ARIMA(0，1，1)(1，1，1)12模型可以較為準(zhǔn)確地預(yù)測(cè)遼寧省AEFI病例報(bào)告趨勢(shì)，可為合理配置AEFI調(diào)查診斷所需資源提供理論依據(jù)。

近年來，隨著擴(kuò)大免疫規(guī)劃工作的深入開展，疫苗可控制傳染病發(fā)病率逐年下降，預(yù)防接種后發(fā)生的疑似預(yù)防接種異常反應(yīng)(adverse event following immunization，AEFI)逐漸進(jìn)入了公眾的視線，成為媒體關(guān)注的重點(diǎn)[1-2]，使得免疫規(guī)劃工作者將越來越多的精力投入到AEFI處置工作中去。為合理科學(xué)地配置AEFI調(diào)查診斷所需資源，需要對(duì)AEFI病例報(bào)告趨勢(shì)進(jìn)行科學(xué)的預(yù)測(cè)。

季節(jié)性自回歸積分滑動(dòng)平均模型(autoregressive integrated moving average model，ARIMA)模型是由Box和Jenkins于20世紀(jì)70年代初提出的一著名時(shí)間序列預(yù)測(cè)方法，所以又稱為Box-Jenkins模型。ARIMA模型的基本思想是：將預(yù)測(cè)對(duì)象隨時(shí)間推移而形成的數(shù)據(jù)序列視為一個(gè)隨機(jī)序列，用一定的數(shù)學(xué)模型來近似描述這個(gè)序列。ARIMA(p，d，q)模型中AR是自回歸，p為自回歸項(xiàng)數(shù)；MA為移動(dòng)平均，q為移動(dòng)平均項(xiàng)數(shù)，d為時(shí)間序列成為平穩(wěn)時(shí)所做的差分次數(shù)。季節(jié)性ARIMA(p，d，q)(P，D，Q)12模型中，P、D、Q分別表示季節(jié)性自回歸、差分和移動(dòng)平均的階次，以大寫與非季節(jié)性的p、d、q區(qū)分[3]。

材料與方法

1.資料來源

以2009-2019年中國(guó)免疫規(guī)劃信息管理系統(tǒng)中遼寧省境內(nèi)報(bào)告AEFI病例為研究對(duì)象，以每月的AEFI報(bào)告數(shù)構(gòu)成時(shí)間序列。

2.季節(jié)性ARIMA模型建模步驟[4]

(1)繪制時(shí)序圖，掌握序列基本趨勢(shì)；

(2)通過純隨機(jī)性檢驗(yàn)考察序列平穩(wěn)性；

(3)通過一階差分和季節(jié)差分對(duì)序列進(jìn)行平穩(wěn)化，通過單位根檢驗(yàn)驗(yàn)證其平穩(wěn)性；

(4)通過觀察平穩(wěn)化后序列的自相關(guān)和偏自相關(guān)圖，對(duì)模型進(jìn)行識(shí)別和優(yōu)化，確定d值、p值、q值和P值，D值，Q值；

(5)通過Ljung-Box檢驗(yàn)等方法檢測(cè)模型殘差，判斷模型的適合性；

(6)利用2019年1-12月AEFI報(bào)告數(shù)據(jù)檢驗(yàn)?zāi)Ｐ偷念A(yù)測(cè)效果。

3.統(tǒng)計(jì)分析軟件

使用R 3.6.1 進(jìn)行數(shù)據(jù)分析。采用utils包中的 read.csv函數(shù)提取原始數(shù)據(jù)，采用stats包中的ts函數(shù)對(duì)原始數(shù)據(jù)進(jìn)行時(shí)間序列處理，用tseries包中的adf.test函數(shù)進(jìn)行單位根檢驗(yàn)，用forecast包中的forecast函數(shù)進(jìn)行預(yù)測(cè)，用stats包中的box.test對(duì)殘差序列進(jìn)行白噪聲檢驗(yàn)(模型有效性檢驗(yàn))。

結(jié) 果

1.時(shí)序圖繪制及檢驗(yàn)

利用stats包中的ts函數(shù)繪制時(shí)序圖(圖1)和時(shí)序分解圖(圖2)，由時(shí)序圖可以看出AEFI人次呈長(zhǎng)期增長(zhǎng)趨勢(shì)，為非平穩(wěn)序列，由時(shí)序分解圖可以看出原數(shù)據(jù)受季節(jié)因素影響。經(jīng)Ljung-Box檢驗(yàn)后，認(rèn)為此序列為非白噪聲序列(P<0.05)。

圖1 2009-2019年遼寧省報(bào)告AEFI人次時(shí)序圖

圖2 2009-2019年遼寧省報(bào)告AEFI人次時(shí)序分解圖

2.時(shí)間序列的平穩(wěn)化

由時(shí)序圖看出原始時(shí)間序列具有上升趨勢(shì)，所以首先通過一階差分消除原序列的趨勢(shì)。隨后為了提取原始時(shí)間序列的季節(jié)性信息，需要對(duì)原始數(shù)據(jù)進(jìn)行一階季節(jié)差分消除季節(jié)趨勢(shì)。經(jīng)單位根檢驗(yàn)，一階季節(jié)差分序列為穩(wěn)態(tài)序列(P<0.05)，可認(rèn)為經(jīng)過一階季節(jié)差分后的時(shí)間序列平穩(wěn)。

3.模型識(shí)別

原序列經(jīng)過一階差分后，序列的趨勢(shì)即消除，得d=1；對(duì)一階季節(jié)差分序列進(jìn)行相關(guān)和偏相關(guān)處理后，得到ACF圖(圖3)和PACF圖(圖4)。觀察圖3，得自相關(guān)系數(shù)一階后未超過±2倍估計(jì)標(biāo)準(zhǔn)差范圍，即自相關(guān)系數(shù)l階以后截尾，初步確定q=1；觀察圖4得偏自相關(guān)系數(shù)3階截尾，初步確定p=3。

圖3 自相關(guān)系數(shù)圖

圖4 偏自相關(guān)系數(shù)圖

4.模型的參數(shù)估計(jì)與檢驗(yàn)

根據(jù)模型參數(shù)檢驗(yàn)結(jié)果和參數(shù)間的相關(guān)系數(shù)，使用auto.arima語句對(duì)模型反復(fù)調(diào)試和檢驗(yàn)，以赤池信息準(zhǔn)則(AIC準(zhǔn)則)作為依據(jù)確定最優(yōu)擬合模型，自動(dòng)輸出相應(yīng)模型階數(shù)，選擇ARIMA(0，1，1)(1，1，1)12模型階數(shù)。對(duì)殘差序列做自相關(guān)圖(圖5)，觀察可知自相關(guān)系數(shù)值在滯后1期迅速衰減，據(jù)此可認(rèn)為殘差序列為白噪聲，信息提取較為充分，符合要求。經(jīng)Ljung-Box檢驗(yàn)后，認(rèn)為此序列為白噪聲序列(P>0.05)，該殘差序列相互獨(dú)立。

圖5 殘差序列的自相關(guān)系數(shù)圖

5.預(yù)測(cè)

利用季節(jié)性ARIMA(0，1，1)(1，1，1)12模型對(duì)遼寧省2019年AEFI報(bào)告數(shù)進(jìn)行預(yù)測(cè)(圖6)，并通過實(shí)際報(bào)告數(shù)對(duì)預(yù)測(cè)結(jié)果進(jìn)行檢驗(yàn)，結(jié)果顯示，該模型的預(yù)測(cè)值均在95%可信區(qū)間內(nèi)，其平均絕對(duì)百分誤差(MAPE)僅在3月時(shí)較高(42.16%)，提示該模型預(yù)測(cè)精度可接受(表1)。

表1 2019年1-12月遼寧省AEFI報(bào)告數(shù)預(yù)測(cè)值與實(shí)際值及檢驗(yàn)

討 論

AEFI處置包括資料的收集、病例的調(diào)查、組織相關(guān)專家的討論、診斷，個(gè)別病例可能需要鑒定，AEFI報(bào)告情況的預(yù)測(cè)[5]，對(duì)AEFI處置資源的合理配置具有重要指導(dǎo)意義。

調(diào)查研究顯示，AEFI的報(bào)告除了疫苗本身及受種者自身情況的影響外，也受到各種不確定因素及措施的影響。季節(jié)性ARIMA時(shí)間序列預(yù)測(cè)模型將各種已知的、未知的因素綜合成統(tǒng)一的影響因素蘊(yùn)含在時(shí)間序列變量中[6]。本文建立季節(jié)性ARIMA模型將遼寧省2009-2019年AEFI報(bào)告數(shù)序列分解為趨勢(shì)、周期、季節(jié)和不規(guī)則四種不可觀測(cè)成份后，將原始時(shí)間序列中隱含的季節(jié)性因素提取出來并予以剔除，消除了不確定因素對(duì)時(shí)間序列的影響。R作為一種開源免費(fèi)的軟件，具有強(qiáng)大的數(shù)據(jù)統(tǒng)計(jì)處理及圖形繪制功能[7]。R軟件中的auto.arima語句，可自動(dòng)根據(jù)AIC準(zhǔn)則，輸出相應(yīng)模型階數(shù)，避免了根據(jù)拖尾和截尾來判斷比較主觀，不夠準(zhǔn)確的問題，提高了模型識(shí)別效率。

本文使用R軟件，經(jīng)序列平穩(wěn)化，季節(jié)趨勢(shì)剔除，參數(shù)估計(jì)、自動(dòng)識(shí)別和檢驗(yàn)等步驟建立了遼寧省AEFI報(bào)告數(shù)季節(jié)性ARIMA模型，通過對(duì)2019年AEFI報(bào)告數(shù)的預(yù)測(cè)及與實(shí)際值的對(duì)比，認(rèn)為該模型可較為有效地?cái)M合AEFI報(bào)告情況，具有一定的推廣價(jià)值。但是，由于AEFI報(bào)告本身具有一定特殊性，如何將其運(yùn)用到實(shí)際中去，仍是一個(gè)值得探討的課題。