ARIMA乘積季節(jié)模型在我國麻疹發(fā)病預(yù)測中的應(yīng)用

羅興甸，戴家佳，羅登菊

(貴州大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院，貴州 貴陽 550025)

麻疹具有很強的傳染性，雖然目前已有麻疹疫苗，但仍有很多幼兒死于麻疹。根據(jù)發(fā)布的相關(guān)報告，2016年全球估計有9萬人死于麻疹，比2000年的55萬例死亡人數(shù)下降84%，標(biāo)志著麻疹死亡數(shù)首次下降到每年10萬例以下。上世紀(jì)中期，我國麻疹年發(fā)病數(shù)曾高達900多萬，到2017年，年發(fā)病人數(shù)僅達6670例，2018年1月至5月，全國麻疹疫情低于近年同期水平。為了實現(xiàn)零麻疹病例，有必要了解我國麻疹發(fā)病規(guī)律，預(yù)測其變化趨勢，據(jù)此作出更好的預(yù)控措施。

許多因素共同影響傳染病的發(fā)病人數(shù)，傳染病的發(fā)生也具有季節(jié)性、趨勢性等特點，所以，靜態(tài)的模型難以揭示其變化趨勢，應(yīng)該根據(jù)發(fā)病數(shù)自身變化規(guī)律建立相應(yīng)的動態(tài)模型。田靈芝等[1]應(yīng)用馬爾可夫模型鏈建立模型分析預(yù)測傳染病數(shù)據(jù)，李若曦等[2]采用ARIMA(autoregressive integrated moving average)模型對北京市豐臺區(qū)流行性感冒的流行性趨勢進行預(yù)測研究，曹磊等[3]利用決策樹模型預(yù)測陜西省手足口病重癥發(fā)病的風(fēng)險及探討可能的危險因素，張喜紅等[4]利用時間序列分析方法中的SARIMA(seasonal autoregressive integrated moving average)模型探討肺結(jié)核的發(fā)病規(guī)律，都取得了較好的效果。

ARIMA和SARIMA技術(shù)已被廣泛應(yīng)用于預(yù)測變量如何隨時間變化。這些技術(shù)通常使用(季節(jié)性)自回歸項和(季節(jié)性)移動平均項來分析及預(yù)測時間序列的變化規(guī)律。如相關(guān)書籍和文獻所記載，這些預(yù)測技術(shù)將變量的先前值和相應(yīng)的誤差項視為預(yù)測未來值的基本信息。運用ARIMA和SARIMA方法能對時間序列數(shù)據(jù)集做出較為準(zhǔn)確的分析預(yù)測。氣溫，電力需求，小麥價格，通貨膨脹，失業(yè)等各個領(lǐng)域的預(yù)測分析證明了ARIMA和SARIMA模型的有效性[5]。本文根據(jù)我國麻疹月發(fā)病數(shù)進行統(tǒng)計建模，文章結(jié)構(gòu)如下：首先介紹了ARIMA乘積季節(jié)模型，然后根據(jù)真實數(shù)據(jù)進行建模，最后進行相關(guān)討論。

1 ARIMA乘積季節(jié)模型

當(dāng)時間序列受到氣候、溫度等隨機效應(yīng)與短期的相關(guān)性交互影響時，隨機效應(yīng)和相關(guān)性存在乘積關(guān)系，通常情況下可以利用ARIMA(p,d,q)(P,D,Q)s乘積季節(jié)模型來探索序列中存在的有用信息。該模型通過d階趨勢差分和D階以周期s為步長的季節(jié)差分運算來使數(shù)據(jù)平穩(wěn)。其表達式為：

(1)

其中，

Φ(Β)=1-φ1Β-…-φpΒp,

Θ(Β)=1-θ1Β-…-θqΒq,

Φs(Β)=1-φ1Βs-…-φPΒPs,

Θs(Β)=1-θ1Βs-…-θQΒQs。

εt表示隨機干擾序列，Φ(Β)dΧt表示同一周期內(nèi)不同周期點的相依關(guān)系，則描述了不同周期的同一周期點上的相關(guān)關(guān)系，二者乘積描述了兩個因素的作用[6]。若P=Q=D=0，則該模型便成為普通的ARIMA模型。

ARIMA模型的建立過程為：首先判斷原始序列是否平穩(wěn)，若不平穩(wěn)，則需要進行預(yù)處理使其平穩(wěn)；然后，將模型中的未知參數(shù)值估計出來；最后對模型進行檢驗[7]。

ARIMA模型建模的前提是平穩(wěn)時間序列，因此，首先利用兩種方法對時間序列進行平穩(wěn)性檢驗，第一種是單位根檢驗，當(dāng)p<0.05時，認為序列平穩(wěn)，反之，則不平穩(wěn)；第二種是圖示法，即通過時序圖檢驗和自相關(guān)圖檢驗作出判斷，若作出的圖呈現(xiàn)序列在一個常值范圍內(nèi)隨機的波動，則說明序列已經(jīng)平穩(wěn)，反之，就不平穩(wěn)。然后根據(jù)平穩(wěn)序列的ACF圖和PACF圖，確定d和D的值，估出模型中p、q、P、Q的值?？筛鶕?jù)所調(diào)查得到的資料確定s的值。對于同一個時間序列，可能會建立多個統(tǒng)計模型，據(jù)AIC準(zhǔn)則選出最佳的擬合模型用于統(tǒng)計推斷。

模型的診斷檢驗主要指檢驗?zāi)Ｐ褪欠裼行?，即：建立的模型是否能充分反映?shù)據(jù)本身包含的信息。若是，則建立的模型有效，若否，則應(yīng)該對模型進行修正調(diào)整使之優(yōu)化。一般用QLB統(tǒng)計量進行檢驗，若p值大于顯著性水平α，則不能拒絕原假設(shè)，即殘差是純隨機序列，說明模型通過檢驗。

最后利用所建立的模型進行短期預(yù)測，進而作出相應(yīng)的評估。

2 模型的建立

2.1 數(shù)據(jù)來源

本文所有數(shù)據(jù)都來自國家衛(wèi)生計生委疾病預(yù)防控制局2004年1月至2018年8月的全部法定報告?zhèn)魅静∫咔橘Y料。將得到的這176個月的我國麻疹月發(fā)病數(shù)構(gòu)成的時間序列記為Χt。

2.2 數(shù)據(jù)的平穩(wěn)性檢驗與處理[8]

采用R軟件進行數(shù)據(jù)的處理和分析[9]。根據(jù)圖示法和單位根檢驗判斷序列是否平穩(wěn)，結(jié)果如下。

圖1 2004年1月至2017年6月我國麻疹月發(fā)病數(shù)(Χt)時間序列圖Fig.1 Time series plot of the monthly number of measles (Χt) from January 2004 to June 2017 in China

圖2 ln(Χt)的時序圖、ACF圖、PACF圖Fig.2 The time series diagram, ACF diagram and PACF diagram of ln(Χt)

由圖1可看出我國麻疹月發(fā)病數(shù)在2004年1月至2017年6月期間以月為周期呈現(xiàn)出規(guī)則的周期性，并且伴隨一定的遞減趨勢，可能因為受到氣候、溫度、濕度等隨機效應(yīng)的影響，每年4月前后全國麻疹發(fā)病數(shù)比其他月份都多。為了消除異方差的影響，對原序列進行對數(shù)變換得到ln(Χt)，圖2的結(jié)果反映序列仍存在明顯的趨勢性，通過單位根檢驗得到p=0.3175>0.05，所以斷定該序列不平穩(wěn)。

為了能更好地用模型對數(shù)據(jù)進行統(tǒng)計分析，對ln(Χt)進行一階趨勢差分和以周期12為步長的季節(jié)差分，獲得新的序列ΔΔ12ln(Χt)，進行單位根檢驗，可以得到p=0.01<0.05，即可認為ΔΔ12ln(Χt)已平穩(wěn)。

2.3 模型的參數(shù)估計與建立

通過以上步驟可確定d=1,D=1,s=12,模型可初步確定為ARIMA(p,1,q)(P,1,Q)12,其中p、q、P和Q是待估的參數(shù)。作出ΔΔ12ln(Χt)的ACF圖與PACF圖如下。

圖3 ΔΔ12ln(Χt)的時序圖、ACF圖、PACF圖Fig.3 The time series diagram, ACF diagram and PACF diagram of ΔΔ12ln(Χt)

從圖3分析可得差分后的麻疹發(fā)病數(shù)呈現(xiàn)隨機波動，除了周期點以外其他階數(shù)的自相關(guān)系數(shù)和偏自相關(guān)系數(shù)都在區(qū)間內(nèi)波動，因此，q=0，但12階函數(shù)不為0，即Q=1；從PACF圖中可得出p=0，P=1或2。因此可供選擇的模型分別為ARIMA(0,1,0)(1,1,1)12和ARIMA(0,1,0)(2,1,1)12。通過R軟件里auto.arima函數(shù)自動選擇階數(shù)進行建模，得到ARIMA(1,0,1)(0,1,1)12。

表1 備選模型參數(shù)估計表Tab.1 Alternative model parameter estimation table

表1反映了三個備選ARIMA模型中的參數(shù)估計值，由表1可得ARIMA(0,1,0)(1,1,1)12模型的AIC等于26.39，ARIMA(0,1,0)(2,1,1)12模型的AIC等于28.28，ARIMA(1,0,1)(0,1,1)12模型的AIC等于25.50，根據(jù)最小信息量準(zhǔn)則(an information criterion,AIC)值最小為最優(yōu)模型選擇標(biāo)準(zhǔn)，選擇ARIMA(1,0,1)(0,1,1)12模型。

首先檢驗?zāi)Ｐ虯RIMA(1,0,1)(0,1,1)12的殘差序列，圖4呈現(xiàn)殘差序列通過了正態(tài)性檢驗；其次，作出模型ARIMA(1,0,1)(0,1,1)12殘差序列的診斷圖(見圖5)，殘差序列的各滯后階數(shù)的自相關(guān)系數(shù)都近似于零，用Ljung-BoxQ統(tǒng)計量進行檢驗，QLB=0.046086，p=0.83>0.05，即模型已全部提取有用信息，所以ARIMA(1,0,1)(0,1,1)12模型通過了殘差檢驗，該模型能很好地擬合麻疹月發(fā)病數(shù)的變化規(guī)律。因此，ARIMA(1,0,1)×(0,1,1)12的表達式為：

(1-0.9667Β)(1-Β12)ln(Χt)

=(1-0.1815Β)(1+0.6614Β12)εt。

(2)

2.4 模型擬合與預(yù)測

首先，利用所建立的ARIMA(1,0,1)×(0,1,1)12模型擬合2004年1月至2017年6月我國麻疹月發(fā)病數(shù)，擬合結(jié)果如圖6所示。圖6展現(xiàn)擬合值與觀測值幾乎重合，表明擬合的效果很好。

圖4 模型ARIMA(1,0,1)×(0,1,1)12的殘差正態(tài)Q-Q圖Fig.4 The residual normal Q-Q plot of the ARIMA(1,0,1)×(0,1,1)12model

圖5 模型ARIMA(1,0,1)×(0,1,1)12的殘差診斷圖Fig.5 The residual diagnosis plot of the ARIMA(1,0,1)×(0,1,1)12model

圖6 ARIMA(1,0,1)×(0,1,1)12模型的擬合結(jié)果圖Fig.6 The fitting result diagram of the ARIMA(1,0,1)×(0,1,1)12model

其次，用ARIMA(1,0,1)×(0,1,1)12模型預(yù)測2017年6月至12月我國麻疹月發(fā)病數(shù)，預(yù)測結(jié)果見表2。由此可得7月之后麻疹月發(fā)病數(shù)的預(yù)測值逐月減少，符合實際情況。相對誤差最大為5.31%，最小為0.23%，分別是在10月和7月，用該模型對麻疹的月發(fā)病數(shù)進行預(yù)測，其相對誤差控制在了5%以內(nèi)，預(yù)測值的平均相對誤差為3.31%，說明取得了較好的短期預(yù)測效果。

最后，將2017年7月至2018年8月的實際數(shù)據(jù)加入原時間序列當(dāng)中重新建立模型，得到ARIMA(1,0,1)×(0,1,2)12模型的表達式為

(1-0.9643Β)(1-Β12)ln(Χt)=(1-0.1812Β)(1+0.606Β12+0.1191Β24)εt。

(3)

用(3)式預(yù)測2018年9月至2019年4月我國的麻疹月發(fā)病數(shù)，得到表3結(jié)果。

3 討論

預(yù)測全國麻疹月發(fā)病情況對如何實施麻疹的預(yù)防與監(jiān)測提供科學(xué)的參考依據(jù)，因此，具有一定的意義。目前國內(nèi)對麻疹的預(yù)測研究都只是針對某一地區(qū)的發(fā)病狀況進行相關(guān)的分析，沒有文獻對全國的麻疹月發(fā)病數(shù)據(jù)進行統(tǒng)計分析。而對于傳染病預(yù)測模型大概包括微分方程模型、余弦模型、灰色預(yù)測模型、Markov模型[10]、人工神經(jīng)網(wǎng)絡(luò)模型[11]、貝葉斯模型等[12]。其中，微分方程模型只描述傳染病在人體內(nèi)的發(fā)展過程，并沒有考慮到人可能在這過程中采取一定的措施和手段去防治傳染病；余弦模型只能對數(shù)據(jù)結(jié)構(gòu)相對簡單的時間序列進行較為準(zhǔn)確的分析；當(dāng)數(shù)據(jù)有較大波動，灰色預(yù)測模型的預(yù)測精確度會大大降低。而ARIMA乘積季節(jié)模型因為能融合序列的短期相關(guān)性和季節(jié)影響因素之間的關(guān)系，從而更能反映數(shù)據(jù)里隱藏的信息，所以廣泛適用于分析較為復(fù)雜的時間序列。

表2 用ARIMA(1,0,1)×(0,1,1)12模型對2017年下半年數(shù)據(jù)預(yù)測結(jié)果Tab.2 The forecast result of using the ARIMA (1,0,1)×(0,1,1)12model to predict the second half of 2017 data

表3 2018年9月至2019年4月我國麻疹月發(fā)病數(shù)預(yù)測值Tab.3 Forecast of the number of measles cases in China from September 2018 to April 2019

從所得數(shù)據(jù)作出的時序圖可以看出，麻疹的月發(fā)病不僅受到季節(jié)的影響，并且隨著時間變化呈現(xiàn)帶有趨勢的變化規(guī)律，所以只考慮一個因素的模型不能對其進行較為準(zhǔn)確的分析預(yù)測。綜合考慮到ARIMA乘積季節(jié)模型的優(yōu)點，通過較為嚴格的統(tǒng)計建模步驟，本文利用該模型對麻疹月發(fā)病數(shù)進行統(tǒng)計分析，還做了短期預(yù)測，取得較為良好的效果。

但是該模型也存在不足的地方：①若獲得的歷史數(shù)據(jù)較少，則所建立的模型有效性相對較低，所以收集得的數(shù)據(jù)越多越能建立更優(yōu)化的模型；②若做短期預(yù)測，效果很好，但做長期預(yù)測，效果不佳。應(yīng)不斷地加入能收集到新的數(shù)據(jù)，以更正或重新擬合更優(yōu)的模型。實際問題的復(fù)雜性會加大建立模型的難度。因此，需要查閱相關(guān)資料全面深入地了解真實情況，研究對數(shù)據(jù)進行預(yù)處理的方法，根據(jù)數(shù)據(jù)本身建立合適的模型，從而做出準(zhǔn)確的分析與預(yù)測。