卞子龍，卓瑩瑩，賀志強，張 楓，蔡奇慧，吳 靜*

1南京醫(yī)科大學生物醫(yī)學工程與信息學院，2公共衛(wèi)生學院，江蘇 南京 211166；3浙江大學公共衛(wèi)生學院，浙江 杭州 310058

肺結(jié)核是由結(jié)核桿菌侵入人體肺部引起的一種慢性呼吸道傳染病，其傳染性強，易反復發(fā)作，可導致癌變惡化的發(fā)生，被列入全球十大致死傳染病之一。據(jù)世界衛(wèi)生組織估計，2017年全球結(jié)核病人數(shù)約為17 億，約有1 000 萬新結(jié)核病患者［1］。中國是全球22個結(jié)核病高負擔國家之一，患者數(shù)量占全球的9%，居全世界第2位［2］。中國國家衛(wèi)生委員會發(fā)布的2017年全國法定傳染病疫情概況顯示，全國（除港澳臺地區(qū)）全年共報告肺結(jié)核發(fā)病835 193例，在乙類傳染病中發(fā)病數(shù)和死亡數(shù)均居第2 位，已經(jīng)成為我國重點關(guān)注的公共衛(wèi)生問題。上海市疾病預防控制中心公布的上海市傳染病疫情報告顯示，2017年全市居民新登記肺結(jié)核3 624 例，發(fā)病率24.9/10萬，較2016 年下降3.0%；外來流動人口新登記肺結(jié)核2 821 例，發(fā)病數(shù)較2016 年下降2.6%。雖然上海肺結(jié)核疫情已得到了有效的控制，但由于耐多藥肺結(jié)核的流行、人口老齡化加速以及外來人口流動性增加等問題，上海的結(jié)核病防治工作又面臨著新的嚴峻考驗，上海地區(qū)的結(jié)核病防控依然不容輕視［3-4］。

時間序列是指將相同統(tǒng)計指標的數(shù)值按其發(fā)生的時間先后順序排列而成的數(shù)列，對其分析的主要目的是根據(jù)對已有歷史數(shù)據(jù)規(guī)律的挖掘從而實現(xiàn)對未來的預測，故近年來被越來越廣泛地應用在傳染病的發(fā)病預測中。本文采用兩種時間序列模型——自回歸移動平均模型（autoregressive integrat?ed moving average model，ARIMA）乘積季節(jié)模型與指數(shù)平滑模型對上海市2007 年1 月—2017 年12 月肺結(jié)核月發(fā)病人數(shù)進行擬合分析，預測2018年1—6月的肺結(jié)核月發(fā)病人數(shù)，并與實際值進行比較，探討這兩種模型在上海市肺結(jié)核疫情預測中的效果，確定預測肺結(jié)核發(fā)病趨勢的最優(yōu)模型，為上海市肺結(jié)核防控工作提供科學依據(jù)。

1 資料和方法

1.1 資料

數(shù)據(jù)資料來源于上海市疾病預防控制中心網(wǎng)站（網(wǎng)址：http：//www.scdc.sh.cn/）2007 年1 月—2018年6 月上海市法定報告?zhèn)魅静∫咔橘Y料，其中2007年1 月—2017 年12 月的肺結(jié)核發(fā)病數(shù)據(jù)用于建立模型，2018年1—6月的數(shù)據(jù)用于驗證模型的預測效果，從而確定最優(yōu)模型。

1.2 方法

1.2.1 ARIMA乘積季節(jié)模型

ARIMA 是由美國統(tǒng)計學家Box 和英國統(tǒng)計學家Jenkins提出的著名時間序列預測模型之一，又稱Box?Jenkins 模型。本研究應用同時帶有季節(jié)性與趨勢性的ARIMA 乘積季節(jié)模型ARIMA（p，d，q）×（P，D，Q）S，其中參數(shù)p、d、q 分別為非季節(jié)自回歸階數(shù)、非季節(jié)差分階數(shù)、非季節(jié)移動平均階數(shù)，P、D、Q分別為季節(jié)自回歸階數(shù)、季節(jié)差分階數(shù)、季節(jié)移動平均階數(shù)，s為季節(jié)周期［5-6］。

ARIMA模型的基本思想是，將預測值隨時間遷移而形成的數(shù)據(jù)序列視為一個隨機序列，用相對應的數(shù)學模型來描述該序列中的自相關(guān)性。當模型被識別后，就可從該時間序列的過去值及現(xiàn)在值來預測未來值。建立ARIMA 時間序列模型可歸納為3個主要步驟：①數(shù)據(jù)的預處理（序列的平穩(wěn)化）：觀察時序圖、自相關(guān)分析圖判斷平穩(wěn)性，通過相應差分進行序列的平穩(wěn)化，進行白噪聲檢驗；②模型的識別、定階與模型參數(shù)估計：采用Box?Jenkins 方法建立ARIMA時間序列分析模型，也就是立足于考察數(shù)據(jù)的樣本自相關(guān)、偏相關(guān)函數(shù)判斷相應的階數(shù)，季節(jié)長度s 可由實際應用背景的分析得到；③模型的診斷檢驗及預測：典型方法是對觀測值和模型擬合值的殘差進行白噪聲分析，同時可以結(jié)合赤池信息準則（Akaike information criterion，AIC）、Schwarz貝葉斯準則（Schwarz Bayesian criterion，SBC），選取較優(yōu)模型進行預測［7］。

1.2.2 指數(shù)平滑模型

指數(shù)平滑法是布朗（Robert G.Brown）提出的一種在移動平均法的基礎上發(fā)展而來的時間序列分析預測方法，通過計算指數(shù)平滑系數(shù)，配合以時間序列預測模型對未來的現(xiàn)象做出預測。事實上，大多數(shù)隨機事件，一般都是近期的數(shù)據(jù)會對現(xiàn)在的影響大，遠期的數(shù)據(jù)會對現(xiàn)在的影響小。指數(shù)平滑法的基本思想就是考慮時間間隔對時間發(fā)展的影響，并且各期權(quán)重隨時間間隔的增大呈指數(shù)衰減。指數(shù)平滑法的預測步驟為：①繪制序列圖；②根據(jù)序列圖確定有效參數(shù)；③繪制擬合曲線圖，并觀察擬合效果；④建立指數(shù)平滑模型，對數(shù)據(jù)進行預測。

根據(jù)序列是否具有長期趨勢與季節(jié)效應，可以把序列分為3 大類，采用3 種不同的指數(shù)平滑模型進行序列預測，具體模型選擇見表1［8］。

表1 指數(shù)平滑預測模型的使用場合Table 1 The usage scenarios of exponential smoothing model

指數(shù)平滑模型含有常規(guī)參數(shù)、趨勢參數(shù)和季節(jié)參數(shù)等3個重要參數(shù)，在通常情況下，應綜合運用整體均值、整體趨勢以及季節(jié)性進行預測，通過不同參數(shù)值的組合進行擬合。在選擇較優(yōu)模型時，通過比較均方根誤差（root mean square error，RMSE）、平均絕對誤差百分比（mean absolute percent error，MAPE）、平均絕對誤差（nean absolute error，MAE）的數(shù)值，綜合選取最優(yōu)模型，并對模型的預測效果進行評價。

1.3 統(tǒng)計學方法

應用SAS 9.4 建立ARIMA 乘積季節(jié)模型；應用R 3.5.0建立Holt?Winters三參數(shù)指數(shù)平滑模型。

2 結(jié)果

2.1 數(shù)據(jù)的初步分析

考慮到上海市經(jīng)濟發(fā)達，人口流動性非常高，因此分析本市居民與外來人口和整體發(fā)病數(shù)之間的關(guān)系。為了直觀地比較，利用軟件畫出本市居民與外來人口隨時間變化的堆積面積圖（圖1）。由圖可知，本市居民與外來人口發(fā)病數(shù)占比隨時間變化相對比較平穩(wěn)，呈季節(jié)性趨勢，也就是整體與部分的病例數(shù)的趨勢相對一致，因此可將上海市的病例數(shù)進行整體分析。

圖1 上海市2007—2017 年本市居民與外來人口肺結(jié)核的發(fā)病數(shù)Figure 1 The number of pulmonary tuberculosis cases among residents and migrants in Shanghai from 2007 to 2017

2.2 ARIMA模型

2.2.1 數(shù)據(jù)預處理

繪制2007 年1 月—2017 年12 月肺結(jié)核發(fā)病數(shù)的時序圖。從圖2 可以觀察到，肺結(jié)核發(fā)病數(shù)隨時間變化總體上呈下降的長期趨勢，并且序列取值以12個月為周期呈現(xiàn)出有規(guī)則的上下波動。具體地，肺結(jié)核發(fā)病數(shù)從每年的1—2月開始上升，在該年的3—4 月先達到1 個高峰，繼而波動式下降，在11—12 月份左右略有上升后再持續(xù)下降到次年1—2 月份，跌落谷底。

圖2 2007—2017年上海市肺結(jié)核月發(fā)病數(shù)時間序列圖Figure 2 Time series of monthly reported number of pul?monary tuberculosis cases in Shanghai from 2007 to 2017

由原始序列圖可知序列不平穩(wěn)，存在周期性。對該序列進行白噪聲檢驗，其自相關(guān)檢查的P值均＜0.05，具有統(tǒng)計學意義，判定上海市肺結(jié)核月發(fā)病數(shù)的時間序列屬于非白噪聲序列。再對序列作1階12 步差分，提取其趨勢效應和季節(jié)效應后，時序圖基本平穩(wěn)（圖3）。

圖3 2007—2017年上海市肺結(jié)核發(fā)病數(shù)差分圖Figure 3 Difference of the number of pulmonary tubercu?losis cases in Shanghai from 2007 to 2017

2.2.2 模型識別與定階

序列具有連續(xù)相關(guān)性和季節(jié)性，說明適合ARI?MA 乘積季節(jié)模型ARIMA（p，d，q）×（P，D，Q）S。經(jīng)1 階12 步差分處理后，序列的長期趨勢和季節(jié)周期性被很好地消除，故判斷d=1，D=1，s=12。根據(jù)差分后序列的自相關(guān)函數(shù)（ACF）圖和偏自相關(guān)函數(shù)（PACF）圖（圖4），ACF 圖顯示延遲1 階自相關(guān)系數(shù)顯著非零，PACF圖顯示延遲1、2階偏自相關(guān)系數(shù)均大于2 倍標準差，故q 可能取0、1，p 可能取0、1、2。此外，考慮序列的季節(jié)自相關(guān)特征，差分后的ACF圖顯示延遲12階自相關(guān)系數(shù)顯著非零，而且延遲24階自相關(guān)系數(shù)也并未完全落入2 倍標準差范圍，由此判斷Q 可能取值為0、1、2。差分后的PACF 圖顯示延遲12階和延遲24階的偏自相關(guān)系數(shù)均顯著非零，故P可能取值為0、1、2。

接著對所有可能合理的模型進行擬合，選出符合建模標準的幾種組合。經(jīng)過幾種模型的比較（表2），選取AIC 為1 375.27，SBC 為1 383.61 均最小的相對最優(yōu)模型ARIMA（0，1，1）×（1，1，2）12，參數(shù)估計均具有統(tǒng)計學意義。

圖4 序列的自相關(guān)函數(shù)（ACF）圖與偏自相關(guān)函數(shù)（PACF）圖Figure 4 Autocorrelation function（ACF）and partial autocorrelation function（PACF）of series

表2 ARIMA（p，d，q）×（P，D，Q）S各種模型比較Table 2 Comparison of various models of ARIMA（p，d，q）×（P，D，Q）S

2.2.3 模型的診斷檢驗及預測

殘差的自相關(guān)檢驗表明序列均為白噪聲，均有P＞0.05，表示應該接受殘差不相關(guān)的零假設，認為殘差序列為純隨機序列，所建立的ARIMA（0，1，1）×（1，1，2）12模型是合適的。此外，從殘差自相關(guān)圖中也能看出殘差序列值都落入置信區(qū)間內(nèi)，符合要求。

2.3 指數(shù)平滑模型

根據(jù)上海市肺結(jié)核月報告發(fā)病例數(shù)序列并結(jié)合指數(shù)平滑法的使用方法，選取Holt?Winters加法指數(shù)平滑法模型、Holt?Winters乘法指數(shù)平滑法模型進行嘗試，對2007 年1 月—2017 年12 月數(shù)據(jù)進行擬合。利用R 3.5.0自動擬合出兩種模型的最優(yōu)參數(shù)，選擇最優(yōu)參數(shù)模型進行比較（表4）。結(jié)果顯示，Holt?Winters 加法指數(shù)平滑法建立模型RMSE、MAPE 均比Holt?Winters乘法指數(shù)平滑法模型的評價參數(shù)小，分別是71.99、14.19%。這表明Holt?Winters 加法指數(shù)平滑法模型比較好。

表3 殘差的自相關(guān)檢查Table 3 Autocorrelation check of residuals

2.4 兩種較優(yōu)模型的對比

采用MAE、RMSE、MAPE 評價兩種較優(yōu)模型的擬合效果（表5）。同時，根據(jù)兩種模型對上海市2018年1—6月肺結(jié)核發(fā)病數(shù)預測值及其相對誤差，比較兩者預測效果（表6）。

由表5 可知Holt?Winters 加法指數(shù)平滑法模型的MAE、MRE、RMSE 均比較小，在本例中該模型要優(yōu)于ARIMA（0，1，1）×（1，1，2）12模型。同時，表6具體展示了兩種模型的預測值與相對誤差，也表明Holt?Winters 加法指數(shù)平滑法模型的總的相對誤差要比ARIMA（0，1，1）×（1，1，2）12模型的小。

最后，繪制這兩種模型的擬合曲線與真實值進行更為直觀的比較（圖5），發(fā)現(xiàn)Holt?Winters 加法指數(shù)平滑模型整體與真實曲線相對擬合效果較好。

表4 最優(yōu)參數(shù)設置和效果評價Table 4 Optimal parameter setting and effect evaluation

表5 兩種模型的擬合效果比較Table 5 Comparison of the fitting effects of the two mod?els

3 討論

肺結(jié)核作為一種成因復雜的慢性傳染病，其發(fā)病率高、病死率高，嚴重威脅人類的健康。隨著人們生活水平的提高與醫(yī)療技術(shù)的進步，結(jié)核病在一定程度上得到了控制，但與此同時，耐多藥肺結(jié)核的流行、人口老齡化加速，又給肺結(jié)核的防治工作帶來了新的挑戰(zhàn)。此外，流動人口對肺結(jié)核防治知識掌握率較低，結(jié)核病就診率低，這也是肺結(jié)核防治工作中的難題［9］。上海市地理位置獨特，經(jīng)濟發(fā)達，人口流動性大，上海地區(qū)結(jié)核病的防治工作尤其需要重視。

表6 兩種模型對預測值的誤差比較Table 6 Comparison of relative errors on predicted values between two models

圖5 實際值與兩種模型擬合值序列比較圖Figure 5 Actual value series and fitted value series of the two models

目前對傳染病疫情進行預測的模型有動態(tài)因子模型、BP人工神經(jīng)網(wǎng)絡模型、馬爾可夫鏈（Markov模型）、廣義回歸神經(jīng)網(wǎng)絡（GRNN）模型、ARIMA 模型等［10-14］。本文所采用的ARIMA 模型和指數(shù)平滑模型都屬于時間序列模型，其基本原理都是通過揭示歷史數(shù)據(jù)隨時間的變化規(guī)律，確定已有時間序列的變化模式，并將這種規(guī)律外延，以此來預測未來現(xiàn)象［12］。ARIMA模型綜合考慮了序列的趨勢、周期的變化及隨機干擾等情況，并用模型參數(shù)進行了量化，能較好地反映時間序列的趨勢和變化，并能判斷季節(jié)效應、趨勢效應和隨機波動等因素。但ARI?MA模型依賴于大量的歷史數(shù)據(jù)，建模過程較復雜，模型的參數(shù)常不容易確定，且要求序列在相當一段時間內(nèi)保持相對的平穩(wěn)才能實現(xiàn)預測的精確有效，故在擬合模型之前，通常需要對非平穩(wěn)序列進行差分等轉(zhuǎn)換。指數(shù)平滑模型是基于對移動平均預測方法的改進，其綜合運用臨近的值、整體趨勢和季節(jié)性來進行預測分析，按均方誤差最小的原則來確定平滑系數(shù)，對近期的值給予更大的權(quán)重，因此近期數(shù)據(jù)對結(jié)果影響較大，而遠期數(shù)據(jù)則影響較小，適用于分析隨時間變化不大的數(shù)據(jù)。相對于ARIMA模型，指數(shù)平滑模型的建模過程更簡單一些［15］。對于上海市肺結(jié)核發(fā)病的歷史數(shù)據(jù)，本文的研究結(jié)果顯示，ARIMA（0，1，1）×（1，1，2）12模型RMSE 為76.27，2018年1—6 月預測值的相對誤差和為0.402；運用指數(shù)平滑法構(gòu)建的最優(yōu)模型是Holt?Winters 加法指數(shù)平滑，其RMSE 為69.61，2018 年1—6 月預測值的相對誤差和為0.292，效果優(yōu)于ARIMA 乘積季節(jié)模型。由此推斷，Holt?Winters加法指數(shù)平滑模型更適合用于上海市肺結(jié)核疫情的短期預測。

ARIMA 乘積季節(jié)模型與指數(shù)平滑模型均以歷史數(shù)據(jù)為基礎，建模前提是數(shù)據(jù)的外延，若外界影響因素突然變化，或是有新變量引入，都會對模型預測效果造成大的影響，降低預測效能。因此這兩個模型更加適用于時間序列的短期預測，對序列的更進一步預測，需要及時更新數(shù)據(jù)資料，添加新的實際值對模型進行修正，然后重新擬合預測。另外，影響肺結(jié)核發(fā)病的因素繁多，欲研究其他因素對肺結(jié)核發(fā)病序列的影響，可以考慮結(jié)合其他模型建立多因素模型共同分析。