基于Python語言的ARIMA模型在天津市結(jié)核病發(fā)病率預測中的應用

張曉卉，姚婷婷，陳 陽，張?zhí)鹛?，馬 駿

(天津醫(yī)科大學公共衛(wèi)生學院流行病與衛(wèi)生統(tǒng)計學系，天津 300041)

結(jié)核病是影響我國人民健康的重大傳染病之一，發(fā)病率和病死率均位于我國傳染病第二位[1-2]。世界衛(wèi)生組織(WHO)2018年報告顯示，結(jié)核病合并免疫系統(tǒng)受損(如感染人類免疫缺陷病毒)、營養(yǎng)不良、糖尿病以及吸煙等將大大提高其病死率[3]。因此，對結(jié)核病發(fā)病率進行預測，根據(jù)其未來預測趨勢為相關(guān)部門制定防控措施，具有實際意義。目前，國內(nèi)已有學者利用時間序列模型對結(jié)核病的發(fā)病趨勢進行預測并取得了較好的效果[4-6]，但對天津市結(jié)核病月發(fā)病率預測時間序列研究鮮有報道。鑒于關(guān)于天津市結(jié)核病發(fā)病率水平的研究較少，本研究基于2004年1月—2016年12月天津市結(jié)核病月發(fā)病率數(shù)據(jù)，通過時間序列分析建立乘積季節(jié)差分自回歸移動平均(autoregressive integrated moving average, ARIMA)模型，擬對天津市結(jié)核病月發(fā)病率進行預測，以期增強結(jié)核病防控預警。當前常見的統(tǒng)計分析軟件為SAS、SPSS、R語言等，三者主要針對統(tǒng)計分析，對于大數(shù)據(jù)的挖掘開發(fā)可視化，用戶普及性等不及Python語言，故本研究應用Python語言進行預測模型的擬合。

1 資料與方法

1.1 資料來源 天津市2004—2016結(jié)核病月發(fā)病率數(shù)據(jù)來源于傳染病網(wǎng)絡直報，數(shù)據(jù)獲取平臺為國家人口與健康科學數(shù)據(jù)共享平臺公共衛(wèi)生科學數(shù)據(jù)中心(http://www.ncmi.cn/info/69/1544)，資料可靠。

1.2 Python語言 Python語言作為當前最熱門的編程語言之一，僅次于java語言和C語言[7]。Python語言不僅可用于統(tǒng)計分析，還被廣泛的應用于數(shù)據(jù)爬取、人工智能等領(lǐng)域[8]，其語言簡單、優(yōu)美，且免費開源。強大的第三方庫支持多種科學計算和統(tǒng)計分析[9]。在時間序列分析中，Python語言建模過程簡單，圖形直觀。當處理的時間序列數(shù)據(jù)量較大時，Python語言可利用其第三方庫pandas，從而規(guī)避循環(huán)，極大地節(jié)省程序運行時間，具有R語言等不具有的自身優(yōu)勢。

1.3 季節(jié)性差分自回歸移動平均(SARIMA)模型 ARIMA模型將數(shù)據(jù)視為隨時間變化的隨機變量，根據(jù)序列的歷史值預測將來值[10]?？紤]到季節(jié)性，使用乘積SARIMA模型。模型常表述為SARIMA(p, d, q) (P, D, Q),其中AR是自回歸，p為自回歸項，MA為移動平均，q為移動平均項數(shù)，d為使時間序列平穩(wěn)的差分次數(shù)，“P,D,Q”為相應帶有季節(jié)性的參數(shù)。

1.3.1 建模步驟[11-14](1)平穩(wěn)性檢驗：依據(jù)時間序列圖、自相關(guān)圖 (autocorrelation function, ACF)、偏自相關(guān)圖 (partial autocorrelation function, PACF)及迪基-福勒檢驗(Dickey-Fuller Test)判斷數(shù)據(jù)是否平穩(wěn)，若不平穩(wěn)需要進行平穩(wěn)性處理。(2)平穩(wěn)性處理：可采取對數(shù)并差分、季節(jié)差分、移動平均等方法使序列平穩(wěn)。(3)白噪聲檢驗：當BOX-Ljung統(tǒng)計量P≤0.05時，判斷序列為非白噪聲序列，可進行后續(xù)分析。(4)定階：根據(jù)差分次數(shù)確定d、D，根據(jù)ACF和PACF確定p、q、P、Q，參考赤池信息準則(Akaike information criterion，AIC)及貝葉斯信息準則(Bayesian information criterion，BIC)確定最優(yōu)模型。(5)參數(shù)估計及殘差分析：確定模型參數(shù)并檢驗其顯著性，對殘差進行白噪聲檢驗。若檢驗通過，則模型建模良好。

1.3.2 模型預測及評價 以天津市2004年1月—2015年12月結(jié)核病月發(fā)病率數(shù)據(jù)作為訓練集擬合模型，以 2016年1—12月數(shù)據(jù)作為驗證集驗證模型擬合精度，預測2017年1月—2019年12月天津市結(jié)核病月發(fā)病率[15]。以誤差均方 (mean square error, MSE)、平均絕對誤差 (mean absolute error, MAE)衡量模型的預測精度和建模效果，其中MSE、MAE越小, 表示預測精度越好, 建模效果越好[16]。

1.4 統(tǒng)計學方法 建模平臺采用基于64位Anaconda(Python3.7)，調(diào)用的模塊主要有numpy、pandas、matplotlib以及statsmodels。

Python建模過程，(1)導入一系列包：import pandas as pd,import numpy as np，import statsmodels.api as sm，import matplotlib.pyplot as plt；(2)導入2004年1月—2015年12月數(shù)據(jù)并將其轉(zhuǎn)換為時間類型數(shù)據(jù)：data=pd.read_csv('jiehebin.csv'),dateparse=lambda dates: pd.datetime.strptime(dates,'%Y-%m'),data=pd.read_csv('jiehebin.csv',parse_dates=['month'],index_col='month',date_parser=dateparse)；(3)差分使序列平穩(wěn)：data_first_difference= data.diff(1)，data_seasonal_first_difference=data_first_difference.diff(12)；(4)模型擬合及預測：Mod=sm.tsa.statespace.SARIMAX(data,trend="n",order=(1,1,1),seasonal_order =(3,1,1,12))，results=mod.fit()，predict_ts=results.predict()。

2 結(jié)果

2.1 流行病學趨勢 以2004年1月—2015年12月數(shù)據(jù)建模，將時間序列圖分解為趨勢性成分、季節(jié)性成分和隨機成分三部分，結(jié)果顯示2004年1月—2015年12月天津市結(jié)核病月發(fā)病率總體呈下降趨勢，2005—2008年出現(xiàn)一個發(fā)病高峰，2009年后大幅度下降，隨后趨于平穩(wěn)。結(jié)核病發(fā)病率于每年3—8月份高發(fā)，冬季驟然降低，提示其具有季節(jié)性。見圖1、2。

圖2 天津市2004年1月—2015年12月結(jié)核病月發(fā)病率時間序列分解圖

2.2 SARIMA模型建模結(jié)果

2.2.1 數(shù)據(jù)預處理 原始時間序列圖經(jīng)Dickey-Fuller Test，結(jié)果顯示P值為0.81，原始序列不平穩(wěn)，需要進行平穩(wěn)性處理。對原始序列進行一次差分和一次季節(jié)差分后序列圖平穩(wěn)，見圖3。ACF和PACF亦顯示數(shù)據(jù)已平穩(wěn)，見圖4。Dickey-Fuller Test結(jié)果顯示P值為0.000029，亦說明序列已平穩(wěn)。經(jīng)差分后的平穩(wěn)序列BOX-Ljung統(tǒng)計量P<0.05，判斷該序列為非白噪聲序列。

圖3 差分后序列時間序列圖及滾動均值、滾動標準差圖

圖4 差分后序列的ACF、PACF

2.2.2 模型識別 由序列一階差分和一階季節(jié)差分初步確定d=1，D=1,初步確定模型形式為SARIMA(p,1,q)×(P,1,Q)12。觀察差分后的ACF和PACF(見圖4)，可見ACF呈截尾，PACF呈拖尾，提示季節(jié)性模型為SARIMA(3，1，1)12模型。觀察SARIMA(1,1,1)×(3,1,1)12殘差序列的ACF和PACF(見圖5)，提示非季節(jié)模型為SARIMA(1，1，1)，故原始序列初步擬合為乘積混合效應模型SARIMA(1,1,1)×(3,1,1)12。為確保篩選出最優(yōu)模型，采用從低階到高階逐個進行嘗試的方法挑選模型參數(shù)[17-18]。初步納入SARIMA(3,1,1)×(3,1,1)12、SARIMA(2,1,1)×(3,1,1)12、SARIMA(1,1,1)×(3,1,1)12、SARIMA(1,1,1)×(3,1,2)12、SARIMA(2,1,1)×(3,1,2)12、SARIMA(3,1,1)×(3,1,2)12六個模型進行試驗，根據(jù)AIC、BIC準則選取其中最小值作為最優(yōu)模型(見表1)，因此原始序列擬合為SARIMA(1,1,1)×(3,1,1)12，擬合后殘差的ACF、PACF見圖5。殘差BOX-Ljung統(tǒng)計量P值為0.493,判斷模型擬合后殘差為白噪聲序列。

表1 擬納入模型AIC、BIC值

圖5 SARIMA(1,1,1)×(3,1,1)12模型擬合后殘差的ACF、PACF

2.2.3 參數(shù)估計及檢驗 模型的參數(shù)估計結(jié)果除ar.L1無統(tǒng)計學意義外，其他參數(shù)均有統(tǒng)計學意義。因此，除去ar.L1，將其他參數(shù)全部列入SARIMA(1,1,1)×(3,1,1)12模型。見表2。

表2 SARIMA(1,1,1)×(3,1,1)12模型的參數(shù)估計

2.2.4 模型擬合 將天津市2004年1月—2015年12月結(jié)核病月發(fā)病率數(shù)據(jù)作為訓練集擬合SARIMA(1,1,1)×(3,1,1)12模型，其中MAE=0.306，MSE=0.224。見圖6。

圖6 2004年1月—2015年12月天津市結(jié)核病月發(fā)病率擬合結(jié)果

2.2.5 模型效果評價 將2016年1—12月結(jié)核病月發(fā)病率進行回代預測，實際發(fā)病例數(shù)均在預測發(fā)病例數(shù)的95%CI內(nèi)，模型擬合良好，具有較好的預測性能，其中MAE=0.169，MSE=0.081，可對天津市結(jié)核病的發(fā)病數(shù)進行較準確的預測。見圖7、表3。

圖7 2016年1—12月天津市結(jié)核病月發(fā)病率預測結(jié)果

表3 2016年1—12月實際發(fā)病率與預測值比較(/10萬)

2.2.6 模型預測 利用SARIMA(1,1,1)×(3,1,1)12對天津市2017年1月—2019年12月肺結(jié)核發(fā)病率進行預測，結(jié)果顯示天津市結(jié)核病月發(fā)病率將總體呈現(xiàn)下降趨勢，春季高發(fā)，冬季發(fā)病率降低，符合結(jié)核病發(fā)病規(guī)律，預測結(jié)果可供參考。見圖8、表4。

圖8 SARIMA對2017年1月—2019年12月天津市結(jié)核病月發(fā)病率預測結(jié)果

表4 SARIMA對2019年7—12月天津市結(jié)核病月發(fā)病的預測值

3 討論

隨著深度學習和數(shù)據(jù)挖掘技術(shù)的日漸成熟，Python語言風靡全球。在信息爆炸和多學科融合的大數(shù)據(jù)時代，Python語言作為一門通用語言在統(tǒng)計學應用中的地位也愈加重要，相較于R語言，Python語言在前期數(shù)據(jù)收集、處理、建模和運行速度等方面顯示出卓越性能，因此本研究使用Python語言建立ARIMA時間序列預測模型。

目前，有多種模型可用于傳染病發(fā)病率的預測，如GM(1,1)模型[19]、馬爾可夫鏈預測模型[20]、ARIMA模型[13,21-22]等。其中ARIMA模型作為最經(jīng)典的預測模型之一，可將時間序列分解為趨勢性成分、季節(jié)性成分和隨機干擾，并對噪聲進行分析處理，預測精度較高，是傳染病時間序列預測模型中最重要的手段[23]。胡曉媛等[24]應用SARIMA模型對全國肺結(jié)核月發(fā)病率進行預測，預測MAE值為0.416992。本研究MAE值為0.169,較之稍高。秘玉清等[25]應用SARIMA模型預測山東省結(jié)核病發(fā)病趨勢，得出發(fā)病率將呈現(xiàn)周期性上升的結(jié)論。鑒于對天津市肺結(jié)核月發(fā)病率研究較少，本研究對天津市結(jié)核病月發(fā)病率的流行病學趨勢作了描述性研究，并預測其未來發(fā)病趨勢，可為天津市相關(guān)部門制定防控措施提供理論依據(jù)。

本研究基于2004年1月—2016年12月天津市結(jié)核病月發(fā)病率資料，將時間序列分解為趨勢性、季節(jié)性和隨機噪聲三部分，分析結(jié)核病的發(fā)病趨勢和季節(jié)性，最終確定SARIMA(1,1,1)×(3,1,1)12為天津市結(jié)核病月發(fā)病率的最終模型。首先，利用 2004年1月—2015年12月數(shù)據(jù)建立最優(yōu)模型，結(jié)果顯示SARIMA(1,1,1)×(3,1,1)12可較準確地擬合實際月發(fā)病率，殘差為白噪聲序列，說明建模良好。然后將2016年1—12月結(jié)核病月發(fā)病率進行回代預測，實際值均在預測值95%CI內(nèi)，說明模型適用于對天津市未來結(jié)核病月發(fā)病率的預測。最后，將SARIMA(1,1,1)×(3,1,1)12模型應用于對2017年1月—2019年12月結(jié)核病月發(fā)病率的預測，可用預測值來估計未來結(jié)核病的流行強度，若2017年1月—2019年12月實際發(fā)病人數(shù)在預測發(fā)病人數(shù)的95%CI內(nèi), 表明當月的結(jié)核病疫情基本正常；若發(fā)現(xiàn)實際發(fā)病人數(shù)處于預測發(fā)病人數(shù)的95%CI外, 則提示結(jié)核病疫情有可能發(fā)生異常[15]。

本研究將2016年1—12月發(fā)病率進行回代預測時，出現(xiàn)8—10月份發(fā)病率相對誤差和其他月份相差較大的情況，其原因：一方面，可能是由于ARIMA模型作為經(jīng)典的線性模型在處理具有非線性特點的時間序列問題上表現(xiàn)出一定的局限性；另一方面，結(jié)核病的發(fā)生發(fā)展具有多因素性，未來考慮將除時間因素以外的其他混雜因素列入模型，以提高其預測精度。

時間序列具有混沌現(xiàn)象，存在內(nèi)在隨機性和不規(guī)則有序性[11]。因此，對時間序列的預測在盡可能抓取線性成分的基礎(chǔ)上還應更多的關(guān)注非線性成分。而ARIMA模型作為一種線性模型在處理非線性成分的問題上具有一定的不足，容易使預測精度降低。目前，關(guān)于基于誤差補償思想和相空間重構(gòu)思想的ARIMA混合模型已經(jīng)嶄露頭角[26-27]。然而，對于ARIMA-SVM、ARIMA-LSTM等混合模型雖預測效果較單純ARIMA模型較好，但解釋性差。在ARIMA模型建模過程中，殘差應通過白噪聲檢驗才能判斷建模是否合理，而對殘差進行相空間重構(gòu)則需要殘差具有混沌性，而非隨機序列，以上矛盾之處此類混合模型無法給出合理的解釋。另外，關(guān)于將ARIMA模型和各類神經(jīng)網(wǎng)絡相結(jié)合的混合模型也層出不窮，研究結(jié)果顯示其預測精度相較于單純ARIMA模型高，但由于神經(jīng)網(wǎng)絡背后數(shù)學理論仍處于“黑箱子”狀態(tài)，其混合模型預測的準確性和重復性仍有待探討。神經(jīng)網(wǎng)絡模型的訓練效果和預測精度隨著數(shù)據(jù)量的增大和數(shù)據(jù)維度的增高而增大，基于時間序列數(shù)據(jù)短期預測的特點，將神經(jīng)網(wǎng)絡運用于時間序列中，其理論有待完善。ARIMA模型常用于短期預測，因此，其預測結(jié)果需要隨著數(shù)據(jù)量的更新而不斷更新，如何將ARIMA和其他各種預測模型如灰色模型、隱馬爾可夫鏈等有效結(jié)合，通過充分提取時間序列的線性和非線性成分，控制隨機誤差，提高預測的準確度和精確度，是今后研究的方向。