基于Python語言的ARIMA模型在天津市結核病發病率預測中的應用

張曉卉，姚婷婷，陳 陽，張甜甜，馬 駿

(天津醫科大學公共衛生學院流行病與衛生統計學系，天津 300041)

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

1 資料與方法

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

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

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

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

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

1.4 統計學方法 建模平臺采用基于64位Anaconda(Python3.7)，調用的模塊主要有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月數據并將其轉換為時間類型數據：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)差分使序列平穩：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 結果

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

圖2 天津市2004年1月—2015年12月結核病月發病率時間序列分解圖

2.2 SARIMA模型建模結果

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

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

圖4 差分后序列的ACF、PACF

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

表1 擬納入模型AIC、BIC值

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

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

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

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

圖6 2004年1月—2015年12月天津市結核病月發病率擬合結果

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

圖7 2016年1—12月天津市結核病月發病率預測結果

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

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

圖8 SARIMA對2017年1月—2019年12月天津市結核病月發病率預測結果

表4 SARIMA對2019年7—12月天津市結核病月發病的預測值

3 討論

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

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

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

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

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