羅興甸，戴家佳，羅登菊

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

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

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

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

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

當時間序列受到氣候、溫度等隨機效應與短期的相關性交互影響時，隨機效應和相關性存在乘積關系，通常情況下可以利用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表示同一周期內不同周期點的相依關系，則描述了不同周期的同一周期點上的相關關系，二者乘積描述了兩個因素的作用[6]。若P=Q=D=0，則該模型便成為普通的ARIMA模型。

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

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

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

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

2 模型的建立

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

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

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

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

圖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ī)則的周期性，并且伴隨一定的遞減趨勢，可能因為受到氣候、溫度、濕度等隨機效應的影響，每年4月前后全國麻疹發(fā)病數(shù)比其他月份都多。為了消除異方差的影響，對原序列進行對數(shù)變換得到ln(Χt)，圖2的結果反映序列仍存在明顯的趨勢性，通過單位根檢驗得到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ù)的自相關系數(shù)和偏自相關系數(shù)都在區(qū)間內波動，因此，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ù)最小信息量準則(an information criterion,AIC)值最小為最優(yōu)模型選擇標準，選擇ARIMA(1,0,1)(0,1,1)12模型。

首先檢驗模型ARIMA(1,0,1)(0,1,1)12的殘差序列，圖4呈現(xiàn)殘差序列通過了正態(tài)性檢驗；其次，作出模型ARIMA(1,0,1)(0,1,1)12殘差序列的診斷圖(見圖5)，殘差序列的各滯后階數(shù)的自相關系數(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 模型擬合與預測

首先，利用所建立的ARIMA(1,0,1)×(0,1,1)12模型擬合2004年1月至2017年6月我國麻疹月發(fā)病數(shù)，擬合結果如圖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模型的擬合結果圖Fig.6 The fitting result diagram of the ARIMA(1,0,1)×(0,1,1)12model

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

最后，將2017年7月至2018年8月的實際數(shù)據(jù)加入原時間序列當中重新建立模型，得到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)式預測2018年9月至2019年4月我國的麻疹月發(fā)病數(shù)，得到表3結果。

3 討論

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

表2 用ARIMA(1,0,1)×(0,1,1)12模型對2017年下半年數(shù)據(jù)預測結果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ù)預測值Tab.3 Forecast of the number of measles cases in China from September 2018 to April 2019

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

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