基于ARIMA和SVM混合模型的流感新增病例預(yù)測

李榮庭，段 鵬，胡 瑞，范蓮靜

(云南民族大學(xué) 數(shù)學(xué)與計(jì)算機(jī)科學(xué)學(xué)院，云南 昆明 650000)

根據(jù)以往的經(jīng)驗(yàn)，我國冬季流感流行時(shí)間一般從12月份開始一直持續(xù)到來年的3—5月.經(jīng)中國疾病預(yù)防控制中心病毒預(yù)防控制所首席專家武桂珍的介紹由于流感和新冠肺炎均可經(jīng)呼吸道飛沫傳播及接觸傳播[1-4]，其均受氣候以及人際接觸等的影響較大：秋冬季節(jié)的低溫天氣會使呼吸道排入外界的病毒存活期延長，且在寒冷刺激下，人的抵抗力相對較弱，病毒趁虛而入，易于形成傳播；人員聚集，則是疫情擴(kuò)大的另一必要條件[5].進(jìn)入秋冬季之后，全國更大范圍的復(fù)工和學(xué)生開學(xué)，會形成廣泛的、不同程度的人員聚集，且秋冬季節(jié)里人們室內(nèi)活動相對頻繁，室內(nèi)通風(fēng)情況不及夏季.加之目前全球范圍內(nèi)許多其他國家的新冠肺炎疫情仍較嚴(yán)峻，尚未得到有效控制等，以上這些因素均可能導(dǎo)致我國在今年冬季面臨新冠肺炎疫情的反彈，其反彈時(shí)間也將很可能與冬季流感流行時(shí)間重疊[6].倘若我們能在流感發(fā)現(xiàn)早期及時(shí)預(yù)測疾病的增長趨勢，很有可能我們就可以采取相應(yīng)的措施，及時(shí)止損.同時(shí)權(quán)威精準(zhǔn)的預(yù)測也可以使人們盡早的把預(yù)防流感重視起來，從而進(jìn)一步減少流感蔓延的規(guī)模，從而也從側(cè)面對新冠疫情的防控起到了積極作用.本課題雖然是以甲型H1N1流感為案例進(jìn)行預(yù)測分析，但搭建構(gòu)造的模型，其實(shí)亦適用與其他類似的傳染病，例如肺結(jié)核、新冠肺炎等.對于后續(xù)的相關(guān)研究，起到了一定的借鑒作用.從另一個(gè)角度出發(fā)，也可以用此模型對未來可能出現(xiàn)的潛在的傳染病威脅進(jìn)行排查.

目前現(xiàn)代醫(yī)學(xué)技術(shù)取得重大發(fā)展，但諸多傳染性疾病仍是人類社會向前發(fā)展的重要阻力之一.查閱全世界范圍內(nèi)已有的流感相關(guān)數(shù)據(jù)進(jìn)行分析研究，多數(shù)預(yù)測模型應(yīng)用多元線性回歸、LASSO回歸以及Ridge回歸模型結(jié)合相關(guān)檢索詞數(shù)據(jù)進(jìn)行建模分析，探討回歸模型與流感疫情預(yù)測的相關(guān)性與可行性[7].常見的傳染病預(yù)測模型還有灰色模型、Markov模型、人工神經(jīng)網(wǎng)絡(luò)模型、通徑分析模型等[8].但這都需要大量數(shù)據(jù)作為支持才能得到理想的結(jié)果，而針對突如其來的大規(guī)模傳染病，其預(yù)測數(shù)據(jù)量小、預(yù)測結(jié)果亟待使用的特點(diǎn)對預(yù)測工作的準(zhǔn)確度和速度提出了雙重要求.基于上述要求，實(shí)驗(yàn)選取少樣本數(shù)據(jù)集，利用數(shù)據(jù)變化的自回歸過程建立模型，使得搭建的模型具有較高的可信度及一定的實(shí)用性.

1 模型介紹

1.1 ARIMA模型

時(shí)間序列是利用統(tǒng)計(jì)學(xué)的基本原理,通過對數(shù)據(jù)的采集選用模型以近似估計(jì),利用模型分析揭示數(shù)據(jù)的內(nèi)在特性[9],達(dá)到推測發(fā)展趨勢規(guī)律的目的[10].自回歸差分移動平均模型(ARIMA)是常用的時(shí)間序列之一.該模型主要分為三個(gè)過程，平穩(wěn)過程(I)、自回歸過程(AR)和移動平均過程(MA).該模型反映了時(shí)間序列過去與現(xiàn)在，未來與現(xiàn)在之間的關(guān)系，適用于短期預(yù)測，需要的數(shù)據(jù)量小、參數(shù)少，這與少樣本數(shù)據(jù)集訓(xùn)練有很高的鍥合度.基于上述特點(diǎn)，ARIMA模型在眾多領(lǐng)域的應(yīng)用都有著不錯(cuò)的變現(xiàn)，尤其是在與時(shí)間相關(guān)的預(yù)測方面[11].

1.2 SVM模型

支持向量機(jī)(support vector machines, SVM)是一種二分類模型，它的基本模型是定義在特征空間上的間隔最大的線性分類器，間隔最大使它有別于感知機(jī)；SVM還包括核技巧，這使它成為實(shí)質(zhì)上的非線性分類器.SVM的的學(xué)習(xí)策略就是間隔最大化，可形式化為一個(gè)求解凸二次規(guī)劃的問題，也等價(jià)于正則化的合頁損失函數(shù)的最小化問題.SVM的的學(xué)習(xí)算法就是求解凸二次規(guī)劃的最優(yōu)化算法[12].相比于傳統(tǒng)的優(yōu)化算法，例如人工蜂群算法、蟻群算法、粒子群算法，SVM的優(yōu)化對象不僅可以是固定參數(shù)也可以是其他數(shù)值.ARIMA模型里的關(guān)鍵參數(shù)都是通過模型自身的差分值進(jìn)行自調(diào)整的，無需再優(yōu)化，對于下文中動態(tài)殘差的優(yōu)化，SVM算法較為合適.

2 基于ARIMA-SVM模型的預(yù)測算法

2.1 ARIMA模型搭建

本文通過將2009年美國甲型流感新增病例數(shù)據(jù)與傳統(tǒng)ARIMA模型結(jié)合[13]，并在此基礎(chǔ)上增加模型評估環(huán)節(jié)以確定階數(shù)，繪制修改后的流程圖如下.

圖1 ARIMA模型構(gòu)建流程示意圖

傳統(tǒng)ARIMA模型在階數(shù)選擇方面，一般采用置信度區(qū)間能夠完全覆蓋的數(shù)值，但是往往沒有容忍度的參數(shù)選擇會使模型錯(cuò)過最佳參數(shù)，故本文加入模型評估，為確定階數(shù)增加一些彈性，使得模型在可接受范圍內(nèi)忽略一些誤差從而找到最佳參數(shù).

2.1.1 平穩(wěn)數(shù)據(jù)

流感的傳播具有聚集性和爆發(fā)性，所以針對數(shù)據(jù)的區(qū)域特點(diǎn)和時(shí)間特點(diǎn)，考慮到模型對數(shù)據(jù)平穩(wěn)性的要求[14]，將美國2009年甲型H1N1新增病例數(shù)數(shù)據(jù)利用差分法進(jìn)行平穩(wěn)，計(jì)算特征點(diǎn)數(shù)據(jù)方差，找到合適的差分階數(shù)，下圖2為差分比較圖.

通過圖2內(nèi)容可以看出，一階差分后數(shù)據(jù)相對穩(wěn)定.又根據(jù)差分后的過濾掉空值缺失數(shù)據(jù)繪制效果展示圖.由圖3可知差分后的平穩(wěn)性要比原數(shù)據(jù)好很多，標(biāo)準(zhǔn)差反映出差分后的數(shù)據(jù)的離散程度在可接受的范圍[15],故確定ARIMA模型的參數(shù)d=1.

圖2 一、二階差分效果展示圖 圖3 差分?jǐn)?shù)據(jù)效果展示圖

2.1.2 自相關(guān)圖和偏自相關(guān)圖的分析

本文利用自相關(guān)函數(shù)(ACF)和偏自相關(guān)函數(shù)(PACF)結(jié)合赤池準(zhǔn)則(AIC)和貝葉斯信息準(zhǔn)則(BIC)進(jìn)行模型評估確定相應(yīng)的階數(shù)p和q，下圖4中的陰影表示置信區(qū)間，可以看出不同階數(shù)自相關(guān)性的變化情況，從而選出p值和q值.p值代表預(yù)測模型中采用的時(shí)序數(shù)據(jù)本身的滯后數(shù)(lags),也叫做AR(auto-regressive)項(xiàng).q值代表預(yù)測模型中采用的預(yù)測誤差的滯后數(shù)(lags)，也叫做MA(moving average)項(xiàng)[9].

針對p值和q值得選取，本模型提出利用Matlab軟件對自相關(guān)函數(shù)圖像和偏自相關(guān)函數(shù)圖像進(jìn)行拉普拉斯(Laplace)濾波處理抽離置信區(qū)間和階數(shù)取值特征點(diǎn)，然后用cftool工具箱中的傅里葉函數(shù)和高斯函數(shù)分別進(jìn)行擬合得到置信區(qū)間邊界函數(shù)s(x)和階數(shù)取值函數(shù)F(x).(由于置信度區(qū)間關(guān)于x軸大致對稱，所以本文選取x軸上方的部分，抽象出置信度區(qū)間邊界函數(shù)s(x)).

圖4 自相關(guān)函數(shù)和偏自相關(guān)函數(shù)選取p、q值

(1)

p值和q值的確定原則是，圖像中置信區(qū)間所涵蓋的最后一個(gè)超出其范圍的滯后數(shù)的取值.但在原有的模型中p值和q值的選取是沒有容忍度的，即只有完全符合置信區(qū)間要求才會被選取.本文在此基礎(chǔ)上提出利用容忍度函數(shù)，如公式(1)，篩選符合條件的最小的p值和q值以保證不會使模型陷入局部最優(yōu)而導(dǎo)致忽略全局最優(yōu).其中xi是待定p或q.根據(jù)上述圖像和公式本確定的p值為3，q值為2.

2.1.3 AR模型與MA模型

根據(jù)之前確定的差分階數(shù)還有p值和q值，建立如下數(shù)學(xué)模型：

(2)

(3)

(4)

其中，yt是當(dāng)前值，μ是常數(shù)項(xiàng)，p是階數(shù)，γi是自相關(guān)系數(shù)，εt是誤差，θi能夠消除波動的最佳參數(shù).從實(shí)際情況出發(fā)由于自回歸模型本身的限制，在數(shù)據(jù)選取過程中既要保證穩(wěn)定性，還要保證自相關(guān)性，如果自相關(guān)系數(shù)φi小于0.5，則不宜采用[13].

2.1.4 模型殘差檢驗(yàn)

基于以上模型，進(jìn)行模型殘差檢驗(yàn)，觀測ARIMA模型方差是否為常數(shù)的正態(tài)分布，以確保模型的可行性和合理性[16].并通過對其準(zhǔn)確的預(yù)測采取下一步相應(yīng)的預(yù)防措施.在選擇了合適的參數(shù)之后，對ARIMA模型的殘差進(jìn)行初步預(yù)測，如圖5.

2.2 SVM模型殘差優(yōu)化

2.2.1 SVM優(yōu)化分析

根據(jù)上文中給出的相關(guān)圖像，不難看出預(yù)測數(shù)據(jù)與實(shí)際數(shù)據(jù)還存在這一定的誤差.故針對預(yù)測中存在的誤差本實(shí)驗(yàn)欲通過支持向量機(jī)(Support Vector Machine,SVM)進(jìn)行優(yōu)化[17].根據(jù)上文，本實(shí)驗(yàn)選取原模型殘差作為對象進(jìn)行優(yōu)化，并綜合兩個(gè)模型的實(shí)驗(yàn)結(jié)果得出混合預(yù)測值作為最終預(yù)測值[11]，圖6為優(yōu)化過程流程圖.對于圖6混合預(yù)測值的詳細(xì)解釋參照下文殘差數(shù)據(jù)導(dǎo)入部分.

圖5 殘差預(yù)測效果圖展示 圖6 加入SVM模型優(yōu)化后的ARIMA模型流程示意圖

具體優(yōu)化的操作步驟主要分為3個(gè)部分，第1個(gè)部分根據(jù)之前ARIMA模型預(yù)測數(shù)據(jù)與實(shí)際數(shù)據(jù)做差得到新的誤差數(shù)據(jù)并利用numpy第3方庫進(jìn)行矩陣化為誤差數(shù)據(jù)貼上時(shí)序標(biāo)簽，如公式(5)；第2部分將得到的新的數(shù)據(jù)再做滑動處理擴(kuò)增數(shù)據(jù)集，滑動次數(shù)按照參照之前ARIMA模型得出的，這樣就得到了新的矩陣，共計(jì)63行5列如圖7所示；第3部分將實(shí)際數(shù)據(jù)的誤差作為標(biāo)簽利用SVM模型進(jìn)行有監(jiān)督的學(xué)習(xí)，具體操作流程參下.

Eerror=yreal-yarima

(5)

2.2.2 導(dǎo)入SVR模塊

由于當(dāng)前環(huán)節(jié)需要處理的數(shù)據(jù)是預(yù)測值與實(shí)際值的誤差，那么自變量就可以定做是出現(xiàn)新增病例之后的每一天，與此同時(shí)因變量也自然而然成了誤差本身.在此基礎(chǔ)上就不難推斷在進(jìn)行優(yōu)化之前需要導(dǎo)入SVR模塊[15]，SVR()就是SVM算法來做回歸用的方法(即輸入標(biāo)簽是連續(xù)值的時(shí)候要用的方法)，通過以下語句來確定SVR的模式(選取比較重要的幾個(gè)參數(shù)進(jìn)行測試).

1) kernel：核函數(shù)的類型

一般常用的有rbf函數(shù)，linear函數(shù)，poly函數(shù)，如圖7所示，將數(shù)據(jù)導(dǎo)入函數(shù)擬合圖像發(fā)現(xiàn)rbf函數(shù)圖像與實(shí)際數(shù)據(jù)擬合度最高，故最佳核函數(shù)為rbf函數(shù)[18].

圖7 誤差數(shù)據(jù)最終形式 圖8 核函數(shù)與數(shù)據(jù)貼合度參考圖

2)C:懲罰因子

C表征對離群點(diǎn)的重視程度，C大小與重視程度成正比，與對誤差分類的懲罰成正比，與margin成反比.當(dāng)C趨近無窮的時(shí)候，表示不允許分類誤差的存在，容易過擬合；當(dāng)C趨于0時(shí)，表示我們不再關(guān)注分類是否正確，容易欠擬合.

3)gamma

gamma是rbf函數(shù)的核系數(shù)且gamma的值必須大于0，其作用是保證模型不會過擬合.

綜上所述SVR模塊的幾個(gè)重要參數(shù)可以由上述規(guī)則基本確定：

SVR(C=10, cache_size=200, coef0=0.0, degree=2, epsilon=0.1, gamma=‘scale’, kernel=‘rbf’, max_iter=-1, shrinking=True, tol=0.001, verbose=False)

2.2.3 殘差值數(shù)據(jù)導(dǎo)入

首先要將每日的預(yù)測值，轉(zhuǎn)化為矩陣形式，方便將數(shù)據(jù)導(dǎo)入SVM.由于目前本實(shí)驗(yàn)無法確定優(yōu)化后的效果優(yōu)秀與否，故要同時(shí)考慮原始ARIMA模型的預(yù)測值yarima與ARIMA-SVM模型的預(yù)測值yarima-svm，如公式4.然后將總預(yù)測值減去實(shí)際值，得到殘差值E，即公式5.在此基礎(chǔ)上將殘差值導(dǎo)入SVR模塊，為了篩選出最合適本數(shù)據(jù)集的C和gamma，還要引入一個(gè)新的變量[12]，記作最小殘差Emin.一方面初始Emin是為了拿來和E進(jìn)行比較，縮小檢索C和gamma的范圍，保證篩選范圍不會太大，導(dǎo)致程序崩潰.另一方面，并將循環(huán)中出現(xiàn)的最小E賦值給Emin，循環(huán)最后找到Emin對應(yīng)的C和gamma，記作Cbest和gammabest流程圖如圖8.

ysum=yarima+yarima-svm，

(6)

(7)

圖9 SVM參數(shù)調(diào)整流程圖

3 實(shí)驗(yàn)與分析

實(shí)驗(yàn)數(shù)據(jù)采用WHO官網(wǎng)公開的美國2009年甲型H1N1流感數(shù)據(jù)集，實(shí)驗(yàn)選擇配置為AMD Ryzen 5中央處理器(主頻 3.4 GHz)，內(nèi)存為8G的電腦為硬件平臺，在 Windows 10，64位操作系統(tǒng)下利用anaconda(python3.7.2)進(jìn)行實(shí)驗(yàn).將上述過程中得到的最佳最佳參數(shù)分別導(dǎo)入到ARIMA和SVM模型中，便可以得到較為準(zhǔn)確的預(yù)測結(jié)果了.

3.1 數(shù)據(jù)集

本文使用的數(shù)據(jù)來自世界衛(wèi)生組織WHO(World Health Organization)公開下載數(shù)據(jù)集.對數(shù)據(jù)進(jìn)行了如下處理：①考慮到地方防疫措施和疫苗推廣對疫情的影響，本實(shí)驗(yàn)選取2009年疫情初期4月24日到7月16日共計(jì) 83 d，美國甲型H1N1新增病例數(shù)信息(其中有 4 d 空值)；②數(shù)據(jù)清洗，取空值前后2天的平均值填充;③將數(shù)據(jù)集按4∶1比例分割為訓(xùn)練集與驗(yàn)證集.

3.2 基于ARIMA和SVM模型的預(yù)測

3.2.1 ARIMA模型預(yù)測結(jié)果

將2009年4月24日以來到7月6日之間美國甲型H1N1流感新增病例數(shù)數(shù)據(jù)導(dǎo)入到已經(jīng)初步搭建好的時(shí)間序列模型，并利用python軟件繪制出模型差分值圖像，如圖9所示.并將差分后的數(shù)據(jù)導(dǎo)入傳統(tǒng)ARIMA模型預(yù)測結(jié)果如圖10.

圖10 ARIMA模型數(shù)據(jù)的差分值 圖11 ARIMA模型預(yù)測效果示意圖

很明顯預(yù)測結(jié)果存在較大偏差，其中原因有可能為以下2點(diǎn)：

1) 在偏離處當(dāng)?shù)叵嚓P(guān)部門采取了有效的防控措施，使得實(shí)際新增病例數(shù)遠(yuǎn)遠(yuǎn)低于預(yù)測值.

2) 模型本身仍存在一定缺陷，需要進(jìn)行進(jìn)一步的優(yōu)化.

基于上述原因，本文提出了利用SVM模型對殘差值進(jìn)行改進(jìn)的方案.

3.2.2 ARIMA-SVM模型優(yōu)化結(jié)果

根據(jù)上文中SVM對殘差值的優(yōu)化，將新的預(yù)測模型預(yù)測值、舊的預(yù)測模型預(yù)測值、真實(shí)的新增病例利用python中的畫圖工具繪制甲流新增病例變化趨勢折線圖12用于進(jìn)行相互比對.

3.2.3 模型結(jié)果評估

為了考察新模型預(yù)測效果，本實(shí)驗(yàn)引入平均絕對百分誤差(MAPE)和均方根誤差(RMSE)2個(gè)指標(biāo)來測定模型預(yù)測的精度.其中MAPE表示預(yù)測結(jié)果每天的預(yù)測結(jié)果與當(dāng)天實(shí)際情況的偏離指數(shù)，RMSE表示總體預(yù)測結(jié)果與實(shí)際情況的偏離指數(shù).故MAPE和 RMSE越小，預(yù)測精度越高.公式與對比見表1：

(8)

(9)

圖12 各模型預(yù)測結(jié)果對比圖

表1 新舊模型殘差評估表

通過上面的表格和對比圖，不難發(fā)現(xiàn)改進(jìn)后的模型精準(zhǔn)度大大增加了，疫情病發(fā)的 50 d 內(nèi)，預(yù)測數(shù)據(jù)與實(shí)際數(shù)據(jù)基本重合，準(zhǔn)確率非常高，50 d 之后預(yù)測與實(shí)際情況開始出現(xiàn)偏差，但走勢相對吻合.

查閱相關(guān)資料2009年6月11日(即疫情出現(xiàn) 50 d 以后)，WHO將甲型H1N1警戒級別從5級提高到最高級6級.這是WHO 40年來第1次把傳染病警戒級別升至最高級別.世衛(wèi)生組織確認(rèn)全球75個(gè)國家和地區(qū)，共確診27737例患者，死亡141例.至此美國人民危機(jī)意識已經(jīng)防疫意識有了本質(zhì)上的覺醒，自覺在家隔離，相關(guān)部門加強(qiáng)衛(wèi)生管理.猜測這可能是導(dǎo)致疫情增長狀況出現(xiàn)拐點(diǎn)與預(yù)測數(shù)據(jù)產(chǎn)生偏差的本質(zhì)原因.正如王菊[19]博士所說，獨(dú)立、透明、有效和及時(shí)的信息流對于國際社會控制不斷出現(xiàn)的疾病至關(guān)重要.隨著甲型HIN1流感的細(xì)節(jié)開始被披露，人們應(yīng)對農(nóng)業(yè)尤其對墨西哥的豬養(yǎng)殖業(yè)進(jìn)行一個(gè)更加全面的評估，這對于為未來吸取經(jīng)驗(yàn)教訓(xùn)至關(guān)重要.當(dāng)然，人們目前的注意力集中在國際社會公共衛(wèi)生反應(yīng)，為疫情作準(zhǔn)備意味著為意外作準(zhǔn)備，在不確定的情況下作出快速和有效的反應(yīng).正如應(yīng)對禽流感的經(jīng)驗(yàn)教訓(xùn)所表明的那樣，這可能需要采取更多的措施，而不只是從上至下的“積極和有力”的技術(shù)應(yīng)對措施.本文的實(shí)證分析為上述衛(wèi)生結(jié)論提供有力證明.

4 結(jié)語

本研究基于2009年在美國和墨西哥等地大面積爆發(fā)的甲型H1N1流感新增病例數(shù)據(jù)，展開關(guān)于傳染性肺炎的研究，通過ARIMA -SVM模型預(yù)測得知甲型H1N1流感新增病例人數(shù)明顯高于實(shí)際值，這可能與當(dāng)年當(dāng)?shù)馗鲊t(yī)療工作者展開的防治工作有關(guān)，本研究結(jié)果也為甲型H1N1流感防治措施的有效性提供了間接證據(jù),為之后的傳染病預(yù)防工作提供參考，從而達(dá)到及時(shí)判斷、抉擇和控制的目的.

根據(jù)4月24日到7月16日數(shù)據(jù)本實(shí)驗(yàn)針對傳統(tǒng)ARIMA模型進(jìn)行改進(jìn)：①提出利用圖像擬合函數(shù)得到了最佳p值和q值；②利用ARIMA預(yù)測值與實(shí)際值之間的殘差序列作為樣本集，導(dǎo)入支持向量機(jī)(SVM)預(yù)測殘差，建立ARIMA-SVM混合模型；③在SVM優(yōu)化過程中加入?yún)?shù)自循環(huán)得到最佳懲罰因子和支持向量數(shù)；④通過比較新的評估參數(shù)MAPE和RMSE，發(fā)現(xiàn)支持向量機(jī)有效擬合了原始序列的非線性部分，ARIMA-SVM混合模型的預(yù)測效果優(yōu)于單一ARIMA模型.

針對本模型可以做的進(jìn)一步研究工作：①優(yōu)化ARIMA模型的數(shù)據(jù)平穩(wěn)過程；②本模型應(yīng)對噪聲的自適應(yīng)性較差，如果在預(yù)測的過程中數(shù)據(jù)被其它人為因素影響會導(dǎo)致預(yù)測誤差變大，故在去噪方面可以進(jìn)一步改進(jìn)；③ARIMA模型研究的對象只與時(shí)間有關(guān)，但在現(xiàn)實(shí)社會中，流行性傳染病的傳播受醫(yī)療水平、文化、經(jīng)濟(jì)等多種因素的影響，在本模型的基礎(chǔ)上可以加入人口演化，醫(yī)療隔離等仿真模型，使得模型精度進(jìn)一步提高.