ARIMA 模型和ARIMA-SVM 模型對(duì)上海市2 型糖尿病患者肺結(jié)核發(fā)病的預(yù)測(cè)效果

黃國(guó)寶 黎衍云 吳 菲 沈 鑫 徐望紅

（1復(fù)旦大學(xué)公共衛(wèi)生學(xué)院流行病學(xué)教研室 上海 200032；2上海市疾病預(yù)防控制中心 上海 200336）

2 型糖尿病（type 2 diabetes mellitus，T2DM）和結(jié)核?。╰uberculosis，TB）是常見的慢性病，兩者高發(fā)年齡吻合，且相互影響。TB 又以肺結(jié)核（pulmonary tuberculosis，PTB）最為常見。流行病學(xué)調(diào)查發(fā)現(xiàn)，T2DM 患者發(fā)生結(jié)核病的概率是非T2DM 患者的2～3 倍，且與地區(qū)、種族無關(guān)［1-2］。我國(guó)首次糖尿病-結(jié)核病雙向篩查結(jié)果顯示，糖尿病患者TB 檢出率為1.25%，遠(yuǎn)高于非糖尿病患者的0.36%。而TB 患者若合并糖尿病，則病程進(jìn)展更快、療效差、耐藥率高、預(yù)后兇險(xiǎn)。上海市嘉定區(qū)2008—2014 年的一項(xiàng)調(diào)查顯示，PTB 合并糖尿病的患者中，痰涂片陽性、合并空洞、治療失敗和轉(zhuǎn)入耐多藥和復(fù)治患者占比分別達(dá)73.67%、56.32%、4.21%和13.16%，遠(yuǎn)高于未合并糖尿病患者的42.51%、52.91%、0.14%和7.02%［3］。鑒于T2DM 患者的TB 發(fā)病率高于一般人群，且兩病共患患者預(yù)后較差，及時(shí)掌握T2DM 人群中TB 發(fā)病情況，將有利于結(jié)核病防控和糖尿病管理。

差分自回歸移動(dòng)平均（autoregressive integrated moving average，ARIMA）模型是預(yù)測(cè)疾病發(fā)病率較為成熟的建模方法，常用于預(yù)測(cè)PTB、HIV、乙肝等傳染病的發(fā)病率［4-5］。中國(guó)是T2DM 與TB 同時(shí)高發(fā)的國(guó)家，T2DM 合并TB 的患者逐年增加?！笆濉比珖?guó)結(jié)核病防治規(guī)劃倡導(dǎo)在糖尿病患者等重點(diǎn)人群中開展TB 的主動(dòng)篩查。因此，急需在數(shù)量龐大的中國(guó)T2DM 患者中建立TB 預(yù)測(cè)模型，以預(yù)測(cè)相關(guān)疾病負(fù)擔(dān)，制定和評(píng)估有針對(duì)性的防控策略和措施。

本研究以上海市2010—2015 年T2DM 患者管理數(shù)據(jù)和結(jié)核發(fā)病監(jiān)測(cè)數(shù)據(jù)為基礎(chǔ)，選擇時(shí)間序列模型中的ARIMA 模型，并用支持向量機(jī)（support vector machine，SVM）對(duì)殘差進(jìn)行擬合，建立結(jié)核病發(fā)病率的組合預(yù)測(cè)模型，對(duì)T2DM 患者的PTB 發(fā)病率進(jìn)行預(yù)測(cè)。

資料和方法

資料來源資料來自2010—2015 年上海市糖尿病管理系統(tǒng)和結(jié)核監(jiān)測(cè)系統(tǒng)數(shù)據(jù)庫。研究對(duì)象納入標(biāo)準(zhǔn)為：（1）上海市戶籍；（2）糖尿病診斷時(shí)間為2010—2015 年；（3）糖尿病確診時(shí)年齡≥35 歲；（4）確診時(shí)未患PTB 或已治愈。共納入T2DM 病例217 999 例。將T2DM 病例的身份證號(hào)與上海市結(jié)核監(jiān)測(cè)系統(tǒng)進(jìn)行鏈接，獲得PTB 發(fā)病信息。在納入上海市糖尿病管理系統(tǒng)和結(jié)核病監(jiān)測(cè)系統(tǒng)時(shí)，所有研究對(duì)象已口頭知情同意所錄入系統(tǒng)信息的分析使用。

ARIMA 模型時(shí)間序列分析的基本思想是將預(yù)測(cè)對(duì)象隨時(shí)間變化而形成的序列看作一個(gè)隨機(jī)序列，除偶然原因出現(xiàn)的極個(gè)別序列值外，時(shí)間序列是依賴于時(shí)間（t）變化的一組隨機(jī)變量，可以利用該序列的過去值、現(xiàn)在值來預(yù)測(cè)序列的未來發(fā)展?fàn)顟B(tài)［6］。ARIMA 模型由Box 和Jenkins 提出，也叫博克斯-詹金斯（簡(jiǎn)稱B-J）法，在ARIMA 模型中預(yù)測(cè)值表達(dá)式為過去若干個(gè)取值和隨機(jī)誤差的線性函數(shù)［7］。其基本結(jié)構(gòu)為不包含季節(jié)性的ARIMA（p，d，q）和包含季節(jié)性的ARIMA（p，d，q）（P，D，Q）。模型建立的步驟如下［8］：

模型識(shí)別 通過時(shí)序圖及單位根檢驗(yàn)判斷序列平穩(wěn)性。非平穩(wěn)的序列一般可以通過差分轉(zhuǎn)換為平穩(wěn)序列，根據(jù)差分的次數(shù)來確定差分（d）和季節(jié)差分（D）。

模型擬合及參數(shù)估計(jì) 根據(jù)序列自相關(guān)函數(shù)（autocorrelation function，ACF）和偏自相關(guān)函數(shù)（partial autocorrelation function，PACF）分析圖［9］，判斷模型自回歸階數(shù)（p）、滑動(dòng)平均階數(shù)（q）、季節(jié)回歸階數(shù)（P）以及季節(jié)滑動(dòng)平均階數(shù)（Q）。

模型檢驗(yàn) 即對(duì)模型的殘差進(jìn)行白噪聲檢驗(yàn)：若白噪聲檢驗(yàn)結(jié)果P＜0.05，則需回到上一步驟，重新進(jìn)行參數(shù)估計(jì)，擬合新的模型。若白噪聲檢驗(yàn)結(jié)果P＞0.05，說明序列信息已經(jīng)提取完整，模型結(jié)果合適，但仍需重新進(jìn)行上一步驟，擬合新的模型。若同時(shí)有幾個(gè)模型殘差序列通過白噪聲檢驗(yàn)，且參數(shù)檢驗(yàn)結(jié)果有意義。則根據(jù)擬合優(yōu)度統(tǒng)計(jì)量標(biāo)準(zhǔn)化貝葉斯判別準(zhǔn)則（standardized bayesian information criteria，SBIC）和赤池信息（Akaike information criterion，AIC）最小的原則選擇最優(yōu)模型。

支持向量機(jī)SVM 是Corinna Cortes 和Vapnik等于1995 年提出的一種機(jī)器學(xué)習(xí)方法，主要用于非線性回歸領(lǐng)域。它能夠更好地解決小樣本集模型難以選擇、經(jīng)驗(yàn)風(fēng)險(xiǎn)原則的樣本容量依賴、維數(shù)災(zāi)難、過擬合而外推差、局部極小點(diǎn)等問題［10］。

ARIMA-SVM 組合預(yù)測(cè)步驟將一組時(shí)間序列的數(shù)據(jù)Yt看成由線性結(jié)構(gòu)Lt和非線性結(jié)構(gòu)Nt組成，即Yt=Lt+Nt。。運(yùn)用ARIMA 模型擬合時(shí)間序列的線性結(jié)構(gòu)，SVM 擬合時(shí)間序列的非線性結(jié)構(gòu)。ARIMA-SVM 具體建模步驟見圖1。

圖1 ARIMA-SVM 的組合預(yù)測(cè)模型流程圖Fig 1 Flow chart of ARIMA-SVM hybrid model establishment

模型的擬合和預(yù)測(cè)效果評(píng)估模型的評(píng)估分為內(nèi)部數(shù)據(jù)評(píng)估（模型擬合程度評(píng)估）和外部數(shù)據(jù)評(píng)估（預(yù)測(cè)效果評(píng)估）。ARIMA 和ARIMA-SVM組合模型的預(yù)測(cè)效果采用近期預(yù)測(cè)法，基于2015 年1 月—12 月PTB 月發(fā)病率數(shù)據(jù)進(jìn)行。為了考察模型預(yù)測(cè)效果，用平均絕對(duì)百分誤差（mean absolute percentage error，MAPE）和均方根誤差（root mean square error，RMSE）來測(cè)定模型預(yù)測(cè)的精度。MAPE 和RMSE 越小，預(yù)測(cè)精度越高。

統(tǒng)計(jì)學(xué)分析運(yùn)用SPSS 23.0 統(tǒng)計(jì)軟件建立和整理發(fā)病資料數(shù)據(jù)庫，利用相關(guān)統(tǒng)計(jì)學(xué)模塊進(jìn)行數(shù)據(jù)處理和分析，實(shí)現(xiàn)ARMIA 建模過程。對(duì)ARIMA模型擬合值與真實(shí)值的殘差序列進(jìn)行樣本重構(gòu)獲得SVM 樣本集，使用R 3.6.2 中的e1071 包對(duì)殘差進(jìn)行SVM 預(yù)測(cè)。tune.svm（）函數(shù)確定核函數(shù)的參數(shù)，使用svm（）函數(shù)擬合模型。

結(jié)果

ARIMA 模型的建立ARIMA 模型建立使用2010 年1 月—2014 年12 月上海市戶籍T2DM 患者PTB 月發(fā)病率，繪制時(shí)間序列圖，并對(duì)其趨勢(shì)性、季節(jié)性和隨機(jī)誤差進(jìn)行分解。由結(jié)果可知，該序列是一組包含季節(jié)和趨勢(shì)的非平穩(wěn)隊(duì)列，初步確定模型結(jié)構(gòu)為ARIMA（p，d，q）（P，D，Q）12（圖2、3）。

圖2 上海市2010—2014 年T2DM 患者PTB 發(fā)病率時(shí)序圖Fig 2 Time series plot of PTB incidence in T2DM patients in Shanghai from 2010 to 2014

圖3 上海市2010—2014 年T2DM 患者PTB發(fā)病率時(shí)間序列分解圖Fig 3 Decomposed time series plot of PTB incidence in T2DM patients in Shanghai from 2010 to 2014

模型識(shí)別對(duì)序列進(jìn)行1 階逐期差分和12 階季節(jié)差分后，時(shí)序圖及單位根檢驗(yàn)（augmented dickey fuller，ADF）結(jié)果顯示，差分后序列呈穩(wěn)定序列（t=-5.789 8，P＜0.01），判斷d=1，D=1。為了判斷模型p、q、P 和Q 值，需要繪制差分處理后的ACF 圖和PACF 圖。由差分后時(shí)間序列的ACF 圖和PACF圖可以看出，PACF 均呈拖尾衰減，ACF 截尾。初步判斷自回歸階數(shù)p=3，滑動(dòng)平均階數(shù)q=0。P、Q的定階比較困難，根據(jù)文獻(xiàn)［11］，模型中的參數(shù)P、Q階數(shù)一般在2 以內(nèi)，故我們選擇0、1、2 進(jìn)行由低階到高階逐個(gè)建模，根據(jù)SBIC 和AIC 最小的原則，選擇最佳模型（圖4、5）。

模型參數(shù)估計(jì)及檢驗(yàn)在嘗試不同模型參數(shù)后，經(jīng)白噪聲檢驗(yàn)和模型參數(shù)檢驗(yàn)獲得3 個(gè)備選模型（表1、2）。

圖4 差分后T2DM 患者PTB 發(fā)病率序列自相關(guān)圖Fig 4 Autocorrelation function graph of differenced PTB incidence series in T2DM patients

圖5 差分后T2DM 患者PTB 發(fā)病率序列偏自相關(guān)圖Fig 5 Partial autocorrelation function graph of differenced PTB incidence series in T2DM patients

表1 備選模型擬合效果及Box-Ljung Q 檢驗(yàn)結(jié)果Tab 1 Fitted effect of alternative models and results of Box-Ljung Q test

表2 備選模型參數(shù)估計(jì)Tab 2 Parameter estimations for alternative models

根據(jù)以上結(jié)果，綜合考慮SBIC 和AIC，最終確定最佳模型為ARIMA（3，1，0）（1，1，0）12（不包含常數(shù)）。殘差序列的Box-Ljung Q 檢驗(yàn)P＞0.05，可以認(rèn)為殘差序列是白噪聲序列，且模型的參數(shù)都有統(tǒng)計(jì)學(xué)意義。模型殘差序列的自相關(guān)系數(shù)和偏自相關(guān)系數(shù)均在95%CI 內(nèi)，進(jìn)一步驗(yàn)證了殘差已是隨機(jī)分布的白噪聲序列，說明所選的模型是恰當(dāng)?shù)模▓D6）。模型表達(dá)式為：

圖6 上海市2010—2014 年T2DM 患者PTB 發(fā)病率殘差序列自相關(guān)圖（左）和偏自相關(guān)圖（右）Fig 6 Autocorrelation（left）and partial autocorrelation（right）graphs of residual sequence of PTB incidence in T2DM patients in Shanghai from 2010 to 2014

ARIMA（3，1，0）（1，1，0）12模型擬合2010—2014年T2DM 患者PTB 發(fā)病率：RMSE=4.05，MAPE=5.8%。

ARIMA-SVM 組合模型將2010—2014 年上海市T2DM 患者PTB 發(fā)病率ARIMA 模型擬合值與實(shí)際值的殘差序列使用SVM 模型進(jìn)行擬合，殘差SVM 預(yù)測(cè)模型參數(shù)：核函數(shù)為radial 函數(shù)，懲罰常數(shù)cost=10，γ=0.1。將SVM 殘差預(yù)測(cè)值加上ARIMA 預(yù)測(cè)值，得到組合模型PTB 發(fā)病率預(yù)測(cè)值。ARIMA 模型和ARIMA-SVM 組合模型得到的2015 年P(guān)TB 發(fā)病率月預(yù)測(cè)值見表3。

表3 2015 年上海市T2DM 患者PTB 月發(fā)病率的兩種模型預(yù)測(cè)值比較Tab 3 Comparison of predicted values of the two models based on monthly incidence of PTB in T2DM patients in Shanghai，2015

模型預(yù)測(cè)效果評(píng)估由表3 可知，ARIMA 模型預(yù)測(cè)上海市2015 年T2DM 患者PTB 發(fā)病率的MAPE 為87%，RMSE 為2.96；ARIMA-SVM 組合模型的MAPE 和RMSE 均小于單純ARIMA 模型，提示ARIMA-SVM 組合模型有效校正了ARIMA模型預(yù)測(cè)值的殘差。但兩個(gè)模型的MAPE 值都較大，預(yù)測(cè)精度都有待提高。

討論

隨著預(yù)測(cè)理論及預(yù)測(cè)技術(shù)的不斷發(fā)展和完善，許多統(tǒng)計(jì)理論、預(yù)測(cè)方法和預(yù)測(cè)模型都被應(yīng)用到傳染病的預(yù)測(cè)中。常見的傳染病預(yù)測(cè)模型有ARIMA模型、灰色模型、Markov 模型、人工神經(jīng)網(wǎng)絡(luò)模型、通徑分析模型等［12］。SVM 也被越來越多地應(yīng)用到農(nóng)業(yè)［13］、通訊［14］、傳染病等領(lǐng)域［15］，與傳統(tǒng)模型相結(jié)合進(jìn)行預(yù)測(cè)，均發(fā)現(xiàn)組合模型的預(yù)測(cè)效果高于單一模型。

本研究基于2010—2014 年上海市戶籍T2DM患者中PTB 發(fā)病人數(shù)，建立并比較了不同參數(shù)的ARIMA 模型預(yù)測(cè)2015 年T2DM 患者PTB 發(fā)病的效果，得到最優(yōu)模型進(jìn)行擬合。利用擬合值與實(shí)際值之間的殘差序列作為樣本集，運(yùn)用SVM 預(yù)測(cè)殘差，建立ARIMA-SVM 混合模型。通過比較兩個(gè)模型的MAPE 和RMSE，發(fā)現(xiàn)SVM 有效擬合了原始序列的非線性部分，ARIMA-SVM 混合模型的預(yù)測(cè)效果優(yōu)于單一ARIMA 模型。

ARMIA 擬合2010—2014 年P(guān)TB 發(fā)病率的效果也較好，比較ARIMA 模型和ARIMA-SVM 混合模型的MAPE 值可以發(fā)現(xiàn)，兩個(gè)模型的預(yù)測(cè)精度還存在較大的提升空間。這可能與PTB 發(fā)病率較低且變化不穩(wěn)定有關(guān)。其次，可能是因?yàn)锳RIMA 分析要求序列是平穩(wěn)序列，雖然可以通過差分將非平穩(wěn)序列轉(zhuǎn)化為平穩(wěn)序列，但是差分意味著信息的損失，導(dǎo)致ARIMA 在預(yù)測(cè)時(shí)誤差較大。此外，時(shí)間序列分析的基本思想是將預(yù)測(cè)對(duì)象隨時(shí)間變化而形成的序列看作一個(gè)隨機(jī)序列，其要求序列的變化只與時(shí)間有關(guān)。在現(xiàn)實(shí)社會(huì)中，PTB 發(fā)病率受經(jīng)濟(jì)、文化、醫(yī)療水平等多種因素的影響。本研究通過ARIMA 模型預(yù)測(cè)得到的PTB 發(fā)病人數(shù)明顯高于實(shí)際值，這可能與近年來PTB 防治工作有關(guān)。本研究結(jié)果可為PTB 防治措施的有效性提供了間接證據(jù)。

上海市糖尿病管理系統(tǒng)和結(jié)核病監(jiān)測(cè)系統(tǒng)建立比較完善，數(shù)據(jù)體量大，資料比較完整，可以進(jìn)一步延長(zhǎng)觀察時(shí)間，利用更長(zhǎng)年份T2DM 患者結(jié)核發(fā)病月數(shù)據(jù)進(jìn)行建模和預(yù)測(cè)，比較不同單一模型和組合模型的擬合效果。需強(qiáng)調(diào)的是，本次預(yù)測(cè)只是從數(shù)據(jù)上反映了傳染病的發(fā)病規(guī)律，有助于預(yù)測(cè)傳染病的發(fā)病趨勢(shì)及負(fù)擔(dān)，但并不能用于探討傳染病發(fā)病的影響因素。當(dāng)前，國(guó)家要求開展重點(diǎn)人群結(jié)核病的篩查，糖尿病患者的結(jié)核病發(fā)病情況備受關(guān)注，一些特異的防控策略和措施可能已經(jīng)影響了該人群結(jié)核病的發(fā)病情況。采用ARIMA-SVM 進(jìn)行預(yù)測(cè)，與觀察到的實(shí)際情況進(jìn)行比較，可在一定程度上反映這些防控策略和措施的效果。

本研究初次比較了單一ARIMA 預(yù)測(cè)模型與組合模型在T2DM 患者PTB 發(fā)病預(yù)測(cè)中的應(yīng)用，研究結(jié)果支持組合模型在該人群中應(yīng)用的有效性。盡管預(yù)測(cè)效果仍然不理想，但在一定程度上提示當(dāng)前PTB 防治措施的效果。將來需要繼續(xù)更新數(shù)據(jù)，優(yōu)化混合模型的參數(shù)，準(zhǔn)確估計(jì)T2DM 患者的PTB發(fā)病情況及疾病負(fù)擔(dān)，指導(dǎo)公共衛(wèi)生工作者做好重點(diǎn)人群的疾病防控工作。