張喜紅，李慧，曹文君，崔永梅

（1. 長(zhǎng)治醫(yī)學(xué)院數(shù)學(xué)教研室，山西 長(zhǎng)治 046000； 2. 北京師范大學(xué)統(tǒng)計(jì)學(xué)院數(shù)理統(tǒng)計(jì)系，北京 100875； 3. 長(zhǎng)治醫(yī)學(xué)院公共衛(wèi)生與預(yù)防醫(yī)學(xué)系流行病與衛(wèi)生統(tǒng)計(jì)學(xué)教研室，山西 長(zhǎng)治 046000； 4. 長(zhǎng)治市疾病預(yù)防控制中心傳染病防控科，山西 長(zhǎng)治 046011）

結(jié)核病是嚴(yán)重危害身體健康的重大傳染病，是世界衛(wèi)生組織和我國(guó)重點(diǎn)控制的傳染病之一。我國(guó)是結(jié)核病高負(fù)擔(dān)國(guó)家，每年新發(fā)結(jié)核病患者約90萬(wàn)例，居全球第3位［1］。最近幾年全國(guó)肺結(jié)核發(fā)病人數(shù)呈持續(xù)下降趨勢(shì)，防治工作效果顯著［1］。季節(jié)自回歸求和滑動(dòng)平均 （seasonal autoregressive integrated moving average，SARIMA） 模型現(xiàn)廣泛應(yīng)用于各種傳染病 （肺結(jié)核、手足口病、流行性腮腺炎、菌痢等［2-5］）的預(yù)測(cè)中。SARIMA模型一般形式為SARIMA（p，d，q）×（P，D，Q）S，其中p是季節(jié)內(nèi)的自回歸階數(shù)，d是季節(jié)內(nèi)的差分運(yùn)算階數(shù)，q是季節(jié)內(nèi)的滑動(dòng)平均數(shù)；P是季節(jié)間自回歸階數(shù)，D是季節(jié)間差分的階數(shù)，Q是季節(jié)間滑動(dòng)平均階數(shù)，s是季節(jié)周期的時(shí)間長(zhǎng)度。長(zhǎng)治市是山西省結(jié)核病高發(fā)區(qū)，為探討肺結(jié)核發(fā)病規(guī)律，本研究對(duì)長(zhǎng)治市 2010年至2017年肺結(jié)核發(fā)病數(shù)進(jìn)行分析，利用時(shí)間序列分析方法中的SARIMA模型探討肺結(jié)核的發(fā)病規(guī)律，為檢測(cè)、預(yù)防和治療肺結(jié)核提供科學(xué)依據(jù)。

1 材料與方法

1.1 資料來(lái)源

通過(guò)“長(zhǎng)治市疾病控制中心傳染病疫情統(tǒng)計(jì)系統(tǒng)”收集長(zhǎng)治市2010年1月至2017年12月肺結(jié)核發(fā)病人數(shù)，利用Excel2007建立原始數(shù)據(jù)庫(kù)，通過(guò)Eviews3.1進(jìn)行統(tǒng)計(jì)分析。

1.2 方法

1.2.1 SARIMA模型簡(jiǎn)介：SARIMA（p，d，q）×（P，D，Q）S也叫乘積SARIMA模型，如果一個(gè)時(shí)間序列既有季節(jié)性又有趨勢(shì)，適合建立的模型為SARIMA模型。構(gòu)建SARIMA模型基本步驟為數(shù)據(jù)預(yù)處理、模型識(shí)別和定階、模型參數(shù)估計(jì)和模型診斷、模型預(yù)測(cè)效果評(píng)價(jià)、預(yù)測(cè)。SARIMA模型的數(shù)學(xué)表達(dá)式為：

φ（B）Φ（BS） （1-B）d（1-BS）DXt=θ（B）Θ（BS） εt

其中

φ（B） =1-φl(shuí)B-…-φpBp

θ（B） =1-θlB-…-θqBq

Φ（BS） =1-ΦlBS-…-ΦpBPS

Θ（BS） =1-ΘlBS-…-ΘQBQS

t表示時(shí)間，Xt表示肺結(jié)核月發(fā)病人數(shù)，B表示滯后算子，εt是白噪聲，φ（B）和Φ （B） 滿足平穩(wěn)性條件，θ （B） 和Θ （B） 滿足可逆性條件，即這4個(gè)多項(xiàng)式的根都在單位圓外。

1.2.2 數(shù)據(jù)預(yù)處理：將2010 年 1 月至 2017 年12月肺結(jié)核月發(fā)病數(shù)數(shù)據(jù)分為2部分：2010 年 1 月至2017 年6月肺結(jié)核月發(fā)病數(shù)數(shù)據(jù)作為樣本，2017年 7 月至 2017 年12月肺結(jié)核月發(fā)病數(shù)數(shù)據(jù)作為預(yù)測(cè)樣本。根據(jù)長(zhǎng)治市 2010 年 1 月至 2017 年 6月肺結(jié)核月發(fā)病數(shù)時(shí)序圖判斷序列變化特征，考察樣本序列是否存在異方差，如果存在異方差做Box-Cox變換，消除異方差；根據(jù)樣本序列自相關(guān)函數(shù) （autocorrelation function，ACF）、偏自相關(guān)函數(shù) （partial autocorrelation function，PACF） 考察樣本序列平穩(wěn)性，經(jīng)過(guò)季節(jié)差分和普通差分后轉(zhuǎn)化為平穩(wěn)序列。

1.2.3 模型的識(shí)別和定階：根據(jù)差分后平穩(wěn)序列ACF 、PACF 函數(shù)的拖尾性或截尾性，對(duì)模型進(jìn)行初步識(shí)別和定階。

1.2.4 模型參數(shù)估計(jì)和模型診斷：根據(jù) AIC （Akaike information criterion） 和SC （Schwarz criterion） 準(zhǔn)則綜合判定模型階數(shù)，利用非線性最小二乘法估計(jì)模型參數(shù)，通過(guò)階數(shù)調(diào)整和參數(shù)顯著性檢驗(yàn)，利用Ljung-BoxQ尾概率 （P值） 在檢驗(yàn)水平α=0.05下，對(duì)模型殘差序列進(jìn)行白噪聲的χ2檢驗(yàn)，判斷建立的SARIMA模型的適合性。

1.2.5 模型預(yù)測(cè)效果評(píng)價(jià)：對(duì)模型的預(yù)測(cè)值進(jìn)行回顧性考核，判定模型的預(yù)測(cè)性能。

1.2.6 預(yù)測(cè)：利用模型預(yù)測(cè)2018年肺結(jié)核發(fā)病情況。

2 結(jié)果

2.1 時(shí)間序列特征分析（圖1）

將2010年1月至2017年6月肺結(jié)核發(fā)病人數(shù)的樣本數(shù)據(jù)記為xt，t=1，2，…，90，曲線圖見(jiàn)圖1A。圖形顯示發(fā)病人數(shù)存在明顯波動(dòng)，1年有1個(gè)高峰，一般在3、4月份；存在異方差、有線性下降趨勢(shì)、具有以年為周期的季節(jié)效應(yīng)。為了消除異方差，對(duì)數(shù)據(jù)做自然對(duì)數(shù)變換［yt=ln （xt） ］，曲線圖見(jiàn)圖1B。圖形顯示，數(shù)據(jù)變換后序列異方差基本消除，線性下降趨勢(shì)明顯，并有周期。

圖1 2010年1月至2017年6月長(zhǎng)治市肺結(jié)核病發(fā)病人數(shù)時(shí)序圖Fig.1 Sequence diagram of the number of patients with pulmonary tuberculosis in Changzhi from January 2010 to June 2017

2.2 ACF分析

時(shí)序圖顯示存在明顯的長(zhǎng)期下降趨勢(shì)和季節(jié)周期性，通過(guò)樣本ACF （圖2） 判斷：ACF較大且為正值，隨著時(shí)滯增加呈現(xiàn)緩慢衰減現(xiàn)象，且在時(shí)滯 12，24等處有峰態(tài)，說(shuō)明該序列具有趨勢(shì)性和季節(jié)性且周期為12 （個(gè)月），故嘗試使用SARIMA模型擬合。

圖2 樣本序列的ACF和PACFFig.2 The ACF and PACF of the sample sequence

2.3 模型識(shí)別和檢驗(yàn)

圖3 新序列的ACFFig.3 The ACF of the new sequence

為了消除趨勢(shì)性，采用 1 階季節(jié)內(nèi)差分對(duì)數(shù)據(jù)yt進(jìn)行處理，得到新序列（1-B） yt樣本ACF函數(shù)拖尾和樣本PACF函數(shù)截尾 （圖3） ，顯示經(jīng)過(guò) 1 階季節(jié)差分處理后的序列季節(jié)內(nèi)、外的趨勢(shì)成分均已消除，但季節(jié)成分仍存在，判斷季節(jié)內(nèi)外的差分階數(shù)分別為d=1，D=0。根據(jù)新序列的ACF和PACF擬定模型的初始階數(shù)范圍，結(jié)合AIC和SC準(zhǔn)則，最終確定沒(méi)有常數(shù)項(xiàng)的SARIMA （2，1，0）×（1，0，1）12為 最 佳 模型，其 AIC = -1.198，SC=-1.086，擬合優(yōu)度R2=0.616，調(diào)整的R2=0.603；得到SARIMA預(yù)測(cè)模型具體表達(dá)式為：

（1-B） （1+0.657B+0.279B2） （1-0.906B12）yt=（1-0.885B12） εt，yt=ln （xt）

其中εt～WN （0，0，1272），見(jiàn)圖4。對(duì)模型殘差進(jìn)行白噪聲檢驗(yàn)，季節(jié)內(nèi)原假設(shè)是殘差序列無(wú)關(guān)，備擇假設(shè)是殘差序列為相關(guān)序列，由檢驗(yàn)結(jié)果得到P ＞0.05，且殘差的ACF和PACF均落在2倍標(biāo)準(zhǔn)差內(nèi)，見(jiàn)圖5。說(shuō)明模型擬合后的殘差序列為白噪聲序列。

圖4 SARIMA 的估計(jì)結(jié)果Fig.4 The estimated results of SARIMA

圖5 殘差序列的ACF和PACFFig.5 The ACF and PACF of the residual sequence

2.4 模型預(yù)測(cè)效果評(píng)價(jià)

應(yīng)用模型擬合長(zhǎng)治市 2010年1月到 2017年6月肺結(jié)核發(fā)病人數(shù)，并對(duì)長(zhǎng)治市 2017 年 7 月至12 月肺結(jié)核發(fā)病人數(shù)預(yù)測(cè)，擬合的肺結(jié)核月發(fā)病人數(shù)和實(shí)際肺結(jié)核月發(fā)病數(shù)具有相同的變化趨勢(shì)，說(shuō)明模型對(duì)該區(qū)的肺結(jié)核月發(fā)病人數(shù)預(yù)測(cè)效果較好。2017年7月至12月實(shí)際月發(fā)數(shù)病人數(shù)與預(yù)測(cè)月發(fā)病人數(shù)一致，模型平均絕對(duì)誤差為6.67，平均相對(duì)誤差為5.96%，預(yù)測(cè)效果較好。見(jiàn)表1。

2.5 預(yù)測(cè)結(jié)果

表1 2017年7月至12月肺結(jié)核發(fā)病人數(shù)實(shí)際值與預(yù)測(cè)值比較Tab.1 Comparison of the actual value with the predictive value of tuberculosis from July 2017 to December 2017

結(jié)果顯示，預(yù)測(cè)2018年總發(fā)病1 326例，1月至12月發(fā)病人數(shù)分別為117、114、127、118、118、108、111、108、100、104、102、100。3月發(fā)病人數(shù)最多，為發(fā)病高峰；1月、4月、6月、7月、9月、10月、12月均比2017年各月發(fā)病人數(shù)多。

3 討論

肺結(jié)核病死亡率高，做好肺結(jié)核防控工作非常重要，肺結(jié)核預(yù)測(cè)是肺結(jié)核防控工作中一個(gè)重要環(huán)節(jié)。文獻(xiàn)資料［6-8］顯示常用預(yù)測(cè)模型有：GM （1，1） 預(yù)測(cè)模型、時(shí)間序列預(yù)測(cè)模型、PB人工神經(jīng)網(wǎng)絡(luò)模型、多元線性回歸、灰色馬爾可夫組合預(yù)測(cè)、ARIMA-GM組合預(yù)測(cè)等。本研究利用長(zhǎng)治市2010年1月到2017年6月肺結(jié)核發(fā)病人數(shù)建立SARIMA預(yù)測(cè)模型，建立的模型為SARIMA 模型，利用模型預(yù)測(cè)2017年7月至12月發(fā)病人數(shù)的平均相對(duì)誤差為5.96%，預(yù)測(cè)效果較好，模型建立較合理。時(shí)間序列預(yù)測(cè)模型方法都有其自身的優(yōu)點(diǎn)和不足，SARIMA 模型是一種短期預(yù)測(cè)精度高的模型。

長(zhǎng)治市是山西省肺結(jié)核病高發(fā)區(qū)，3月份為發(fā)病高峰時(shí)期，近幾年年發(fā)病人數(shù)呈明顯下降趨勢(shì)，但2017年發(fā)病人數(shù)比2016年略有增加，預(yù)測(cè)結(jié)果顯示2018年發(fā)病人數(shù)也比2017年稍多。

綜上所述，本研究建立了時(shí)間序列模型SARIMA 來(lái)總結(jié)長(zhǎng)治市肺結(jié)核的發(fā)病規(guī)律，并有效預(yù)測(cè)肺結(jié)核發(fā)病人數(shù)。時(shí)間序列SARIMA預(yù)測(cè)模型適合短期預(yù)測(cè)，長(zhǎng)期預(yù)測(cè)不僅需要調(diào)整數(shù)據(jù)和模型參數(shù)，也要考慮一些突發(fā)素、環(huán)境因素等，這樣才能使預(yù)測(cè)結(jié)果更加精準(zhǔn)。

［1］ 國(guó)務(wù)院辦公廳關(guān)于印發(fā)”十三五”全國(guó)結(jié)核病防治規(guī)劃的通知［EB/OL］.［2017-02-16］.http：//www.dzwww.com/xinwen/guoneixinwen/201702/t20170216_15551354.htm.

［2］ 孟蕾，王玉明.ARIMA模型在肺結(jié)核發(fā)病預(yù)測(cè)中的應(yīng)用［J］.中國(guó)衛(wèi)生統(tǒng)計(jì)，2010，27 （5）：507-509，DOI：10.3969/j.issn.1002-3674.2010.05.018.

［3］ 楊仁東，胡世雄，鄧志紅，等. 湖南省手足口病發(fā)病趨SARIMA模型預(yù)測(cè)［J］. 中國(guó)公共衛(wèi)生，2016，32 （1）：48-52.DOI：10.11847/zgggws2016-32-01-15.

［4］ 李潤(rùn)滋，章濤，梁玉民，等. SARIMA 模型在流行性腮腺炎發(fā)病預(yù)測(cè)中的應(yīng)用［J］. 山東大學(xué)學(xué)報(bào) （醫(yī)學(xué)版） ，2016，54 （9）：82-86.DOI：10.6040/j.issn.1671-7554.0.2015.1163.

［5］ 胡建利，梁祁，吳瑩，等. 季節(jié)時(shí)間序列模型在菌痢發(fā)病預(yù)測(cè)中的應(yīng)用［J］. 中國(guó)衛(wèi)生統(tǒng)計(jì)，2012，29 （1）：34-39.DOI：10.3969/j.issn.1002-3674.2012.01.010.

［6］ 柳巍，曾令城，李煥芝，等. 灰色殘差 GM （1，1） 模型在預(yù)測(cè)肺結(jié)核流行趨勢(shì)中的應(yīng)用［J］. 河南醫(yī)學(xué)研究，2015，24 （7）：1-3.DOI：10.3969/j.issn.1004-437X.2015.07.001.

［7］ 金如鋒，黃成鋼，邱宏，等. 4種模型對(duì)我國(guó)某地區(qū)肺結(jié)核發(fā)病率的預(yù)測(cè)［J］. 現(xiàn)代預(yù)防醫(yī)學(xué)，2008，35 （4）：4866-4869.DOI：10.3969/j.issn.1003-8507.2008.24.047.

［8］ 易靜，胡代玉，楊德香，等. 三種預(yù)測(cè)模型在肺結(jié)核發(fā)病預(yù)測(cè)中的應(yīng)用［J］. 中國(guó)全科醫(yī)學(xué)雜志，2012，15 （5A）：1495-1497.DOI：10.3969/j.issn.1007-9572.2012.13.022.