黃國寶 黎衍云 吳 菲 沈 鑫 徐望紅

（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］。我國首次糖尿病-結(jié)核病雙向篩查結(jié)果顯示，糖尿病患者TB 檢出率為1.25%，遠(yuǎn)高于非糖尿病患者的0.36%。而TB 患者若合并糖尿病，則病程進展更快、療效差、耐藥率高、預(yù)后兇險。上海市嘉定區(qū)2008—2014 年的一項調(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ù)后較差，及時掌握T2DM 人群中TB 發(fā)病情況，將有利于結(jié)核病防控和糖尿病管理。

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

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

資料和方法

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

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

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

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

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

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

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

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

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

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

結(jié)果

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

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

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

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

模型參數(shù)估計及檢驗在嘗試不同模型參數(shù)后，經(jīng)白噪聲檢驗和模型參數(shù)檢驗獲得3 個備選模型（表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 檢驗結(jié)果Tab 1 Fitted effect of alternative models and results of Box-Ljung Q test

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

根據(jù)以上結(jié)果，綜合考慮SBIC 和AIC，最終確定最佳模型為ARIMA（3，1，0）（1，1，0）12（不包含常數(shù)）。殘差序列的Box-Ljung Q 檢驗P＞0.05，可以認(rèn)為殘差序列是白噪聲序列，且模型的參數(shù)都有統(tǒng)計學(xué)意義。模型殘差序列的自相關(guān)系數(shù)和偏自相關(guān)系數(shù)均在95%CI 內(nèi)，進一步驗證了殘差已是隨機分布的白噪聲序列，說明所選的模型是恰當(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 模型擬合值與實際值的殘差序列使用SVM 模型進行擬合，殘差SVM 預(yù)測模型參數(shù)：核函數(shù)為radial 函數(shù)，懲罰常數(shù)cost=10，γ=0.1。將SVM 殘差預(yù)測值加上ARIMA 預(yù)測值，得到組合模型PTB 發(fā)病率預(yù)測值。ARIMA 模型和ARIMA-SVM 組合模型得到的2015 年P(guān)TB 發(fā)病率月預(yù)測值見表3。

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

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

討論

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

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

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

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

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