陳 寶，丘美嬌，林尤斌，周小莉，吳 彪，b

(海南省人民醫(yī)院，海南醫(yī)學(xué)院附屬海南醫(yī)院a.疾病預(yù)防控制辦公室； b.感染科，?？?570311)

流行性感冒(以下簡稱流感)是一種急性呼吸道傳染病[1]，其病原體為病毒，該病毒通過呼吸道傳播，由于其具有較強(qiáng)的傳染性，特別是甲型流感病毒已發(fā)生變異，人群普遍易感，發(fā)病率高，在世界范圍內(nèi)引起多次暴發(fā)流行。在WHO的推動下，流感成為第一個全球監(jiān)測的傳染性疾病，我國的流感監(jiān)測工作始于1952年，1957年成立國家流感中心[2]，我國的流感中心也是全球流感監(jiān)測與應(yīng)對網(wǎng)絡(luò)(Global Influenza Surveillance and Response System，GISRS)的一員[3]，組織和管理著全國的流感監(jiān)測系統(tǒng)，開展流感樣病例(influenza-like illness，ILI)監(jiān)測和流感樣病例暴發(fā)疫情監(jiān)測工作[4]。本研究通過海南省國家級流感監(jiān)測哨點(diǎn)醫(yī)院2016年第1周至2019年第52周的流感樣病例監(jiān)測數(shù)據(jù)，建立預(yù)測流感樣病例發(fā)病趨勢的最優(yōu)季節(jié)性自回歸滑動平均混合(seasonal autoregressive integrated moving average，SARIMA)模型，為當(dāng)?shù)乜茖W(xué)防治ILI提供依據(jù)。

1 模型擬合的基本思想

SARIMA模型是在傳統(tǒng)自回歸滑動平均混合(ARIMA)模型上的推廣，SARIMA模型的簡略寫法是運(yùn)用“后移算子”(B)，B意味著把序列后移一時間點(diǎn)進(jìn)行觀察，若時間序列為Z，則B(Zt)表示Zt-1，當(dāng)進(jìn)行季節(jié)性后移時，在B的右上角注明所需的滑動次數(shù)即可，如B12(Z)為Z的12個觀察之前的值——即季節(jié)性后移變換。

SARIMA模型記為SARIMA(p，d，q)(P，D，Q)s，P、D、Q、s分別為季節(jié)周期自回歸的階數(shù)、差分次數(shù)、滑動平均的階數(shù)、循環(huán)的周期數(shù)。模型的參數(shù)d值、D值為使序列達(dá)到平穩(wěn)時的差分次數(shù)；根據(jù)序列平穩(wěn)時的自相關(guān)函數(shù)(autocorrelation function，ACF)分析圖、偏自相關(guān)函數(shù)(partial autocorrelation function，PACF)分析圖確定p值和q值；再根據(jù)殘差自相關(guān)函數(shù)分析圖和殘差偏自相關(guān)函數(shù)分析圖確定P值和Q值[5-6]。

SARIMA模型的一般形式為：

Φ(BS)(φB)(1-B)d(1-BS)DZt=Θ(BS)θ(B)εt，φ(B)=1-φ1B-φ2B2-…-φpBp

θ(B)=1-θ1B-θ2B2-…-θqBq

Φ(Bs)=1-φs，1Bs-φs，2B2S-…-φs，PBPs

Θ(Bs)=1-φs，1BS-φs，2B2s-…-φs，QBQs

2 資料與方法

2.1 相關(guān)定義

ILI：發(fā)熱(體溫≥38.0 ℃)，伴咳嗽或咽痛之一者[7]。

ILI百分比(ILI%)，為ILI病例占門急診病例的百分比：ILI病例總數(shù)/門急診病例總數(shù)×100%。

2.2 資料來源

本研究采用的研究資料來源于海南省國家級流感監(jiān)測哨點(diǎn)醫(yī)院(海南省人民醫(yī)院)每周一報告的2016年第1周至2019年第52周的流感樣病例數(shù)。海南省人民醫(yī)院全年開展ILI監(jiān)測工作，監(jiān)測科室每日登記匯總流感樣病例數(shù)，每周一24時前經(jīng)醫(yī)院疾控辦審核原始登記資料后通過中國流感監(jiān)測信息系統(tǒng)上報。省、市、區(qū)三級疾控部門定期對哨點(diǎn)醫(yī)院開展督導(dǎo)工作，確保數(shù)據(jù)的真實(shí)性、完整性、準(zhǔn)確性。

2.3 統(tǒng)計學(xué)方法

利用2016年第1周至2019年第52周每周的ILI%數(shù)據(jù)擬合自回歸滑動平均混合(ARIMA)模型。應(yīng)用擬合優(yōu)度R2值、BIC值、Ljung-Box檢驗(yàn)對預(yù)測模型進(jìn)行評估，以確定最優(yōu)模型。采用游程檢驗(yàn)判斷時間序列的平穩(wěn)性，檢驗(yàn)水準(zhǔn)取α=0.05。所有模型分析和擬合均通過SPSS 25.0軟件進(jìn)行。

3 結(jié)果

3.1 2016年第1周至2019年第52周的ILI%情況

2016年第1周至2019年第52周的ILI%數(shù)據(jù)時間序列圖顯示，該院ILI%的原始序列呈現(xiàn)“雙峰”特點(diǎn)，存在季節(jié)波動，主要集中在每年的2月和10月，該院ILI%的原始時間序列為非平穩(wěn)序列(Z=-9.366，P<0.001)，存在明顯的波動趨勢，見圖1。

時間(年-周)圖1 2016年第1周至2019年第52周的ILI%原始時間序列圖

3.2 ARIMA模型擬合

3.2.1 序列的平穩(wěn)化

由圖1可知，該院ILI%的原始時間序列為非平穩(wěn)序列，以1年(52周)為一周期，具有一定的季節(jié)性，經(jīng)1次差分、1階周期為52的季節(jié)性差分，新的ILI%時間序列均值在0刻度線上下波動，基本認(rèn)為該時間序列平穩(wěn)化(圖2)，即周期s=52，d=1，D=1，模型可初步定為ARIMA(p，1，q)(P，1，Q)52。

時間(年-周)轉(zhuǎn)換:自然對數(shù),差異(1),季節(jié)性差異(1,周期52)。圖2 平穩(wěn)后的ILI%序列圖

3.2.2 模型的識別與定階

平穩(wěn)后新序列的ACF圖、PACF圖分別顯示自相關(guān)系數(shù)在1階截尾，可確定p=1、q=1。則模型進(jìn)一步定為ARIMA(1，1，1)(P，1，Q)52。通常情況下，P、Q值不會超過2，且不會同時為0[8]，初步判斷模型為：模型1、ARIMA(1，1，1)(0，1，1)52；模型2、ARIMA(1，1，1)(1，1，0)52；模型3、ARIMA(1，1，1)(1，1，1)52。各備選模型的擬合優(yōu)度檢驗(yàn)結(jié)果(表1)顯示，模型1ARIMA(1，1，1)(0，1，1)52的R2值最接近1、BIC值最小，且Ljung-Box檢驗(yàn)Q值為14.690，P>0.05，即顯示該模型殘差屬于白噪聲，且模型1的殘差序列ACF和PACF均在可信區(qū)間內(nèi)(圖3)，故最合適的模型為ARIMA(1，1，1)(0，1，1)52。

表1 3個備選ARIMA模型的擬合優(yōu)度檢驗(yàn)和假設(shè)檢驗(yàn)

圖3 ARIMA(1，1，1)(0，1，1)52模型殘差序列的自相關(guān)和偏自相關(guān)圖

3.2.3 模型擬合與預(yù)測

應(yīng)用ARIMA(1，1，1)(0，1，1)52模型獲得該院2016年第1周至2019年第52周ILI%的實(shí)測值與擬合值基本吻合(圖4)。

時間(年-周)圖4 ARIMA(1，1，1)(0，1，1)52模型實(shí)測值與擬合值的比較

使用ARIMA(1，1，1)(0，1，1)52模型預(yù)測2019年第41—52周的ILI%，與實(shí)際值對比發(fā)現(xiàn)，該模型在第41-44周短期內(nèi)的預(yù)測的相對誤差較小，而長期預(yù)測的相對誤差較大，見表2。

表2 2019年第41—52周實(shí)際ILI%與預(yù)測值比較

4 討論

目前，越來越多的數(shù)學(xué)模型被用于傳染病預(yù)測，如灰色模型[9]、神經(jīng)網(wǎng)絡(luò)模型[10]、ARIMA模型[11]、SARIMA模型[10]、傳染病動力學(xué)模型[12]。不同的模型各有其特點(diǎn)及使用范圍。在本研究中，該院ILI在每年的2月及10月出現(xiàn)高峰，與國內(nèi)研究基本一致[13-14]。ILI受季節(jié)影響較明顯，故本研究采用季節(jié)性ARIMA模型。

通過分析流感樣病例監(jiān)測指標(biāo)，能夠提早1～2周預(yù)測流感疫情[15]，因此對流感的預(yù)防和控制具有重要的意義。本研究通過2016年第1周至2019年第52周海南省國家級流感監(jiān)測哨點(diǎn)醫(yī)院ILI監(jiān)測數(shù)據(jù)建立季節(jié)性ARIMA模型，經(jīng)序列平穩(wěn)后等一系列處理后，選取ARIMA(1，1，1)(0，1，1)52為最優(yōu)模型，且模型擬合結(jié)果顯示實(shí)測值與預(yù)測值基本吻合，表明該模型對該院ILI%的短期流行趨勢具有較高的預(yù)測精度。制定ILI%預(yù)測就診率及95%可信區(qū)間范圍，有利于醫(yī)院疾病防治部門提升ILI%監(jiān)測工作質(zhì)量及優(yōu)化臨床診治工作：1)當(dāng)實(shí)際發(fā)病率在預(yù)測值的95%可信區(qū)間范圍內(nèi)，則表明當(dāng)下疫情處于一般水平；2)當(dāng)實(shí)際發(fā)病率超過預(yù)測值的95%可信區(qū)間上限時，則應(yīng)警惕流感爆發(fā)或流行，應(yīng)提早部署發(fā)熱門診等醫(yī)護(hù)人員，妥善開展ILI病例的診治工作；3)當(dāng)實(shí)際發(fā)病率低于預(yù)測值的95%可信區(qū)間下限時，則提示可能存在病例信息漏報等行為。

ARIMA模型只需要發(fā)病率即可進(jìn)行預(yù)測，本研究通過ILI%建立ARIMA模型發(fā)現(xiàn)，其短期預(yù)測的相對誤差較小，而長期預(yù)測的相對誤差較大，這是由于時間序列這種精度較高的短期預(yù)測模型隨著時間推移，其預(yù)測精度會降低。這也從側(cè)面反映出ARIMA模型自身方面的不足，原因在于其基于數(shù)據(jù)本身的一種方法學(xué)，沒有考慮到諸如行政手段干預(yù)等一系列社會因素的影響，即只能進(jìn)行短期預(yù)測。但ARIMA模型通過整合趨勢、周期和隨機(jī)誤差等因素，包含了回歸分析和移動平均的優(yōu)點(diǎn)，仍舊成為時間序列短期預(yù)測最為成熟的方法之一。

盡管ARIMA模型構(gòu)建過程中可參照ACF圖及PACF圖選取模型指標(biāo)，但其模型構(gòu)建很大程度上也和研究者自身的經(jīng)驗(yàn)相關(guān)，研究者所構(gòu)建的所謂“最優(yōu)模型”實(shí)質(zhì)上是基于自身經(jīng)驗(yàn)得到的“最優(yōu)選擇”[16]，盡管相關(guān)領(lǐng)域的專家指出p、q值一般不會超過5，但有研究顯示p值為7時，其模型的絕對誤差和相關(guān)誤差更小，擬合效果更好[17]。

綜上，在今后的研究中，盡可能全面收集影響ILI發(fā)病的因素，如氣象因素、經(jīng)濟(jì)水平、地理位置、病毒變異等，實(shí)時動態(tài)調(diào)整模型數(shù)據(jù)，提高模型精度，選取最優(yōu)的ARIMA模型進(jìn)行流感樣病例監(jiān)測及預(yù)測，進(jìn)一步借助信息化系統(tǒng)，及時發(fā)布ILI發(fā)病高峰預(yù)警，高效開展ILI防治工作。