張世潔，溫 瑩，祝 方，程錦泉1，

1)鄭州大學(xué)公共衛(wèi)生學(xué)院流行病學(xué)教研室 鄭州 450001 2)深圳市疾病預(yù)防控制中心 廣東深圳 518073

急性出血性結(jié)膜炎(acute hemorrhagic conjunctivitis,AHC)俗稱紅眼病，是《中華人民共和國(guó)傳染病防治法》中規(guī)定報(bào)告的一種丙類傳染病。AHC是一種具有高度傳染性的病毒性疾病，由腸道病毒70(EV70)、柯薩奇病毒A24變種(CA24v)或腺病毒引起[1-2]。該病潛伏期短，起病急，臨床表現(xiàn)主要為雙眼疼痛、眼瞼腫脹、結(jié)膜充血、角膜炎、異物感和眼部分泌物增多[3]，其傳染性強(qiáng)，主要通過(guò)水或直接接觸傳播，人群普遍易感，較易出現(xiàn)暴發(fā)流行[4]。該病一年四季均可發(fā)生，具有周期性和季節(jié)性特點(diǎn)，夏秋季高溫濕熱的氣候環(huán)境有利于病原體的繁殖與傳播，多為發(fā)病的高峰期[5]。先前關(guān)于AHC的多數(shù)研究都集中在流行病學(xué)和病因?qū)W特征上，少有AHC的時(shí)間序列分析研究。本文利用Python軟件建立季節(jié)性差分自回歸求和滑動(dòng)平均(seasonal auto-regressive integrated moving average,SARIMA)模型，對(duì)河南省AHC的發(fā)病進(jìn)行分析和預(yù)測(cè)，以期了解AHC短期發(fā)病趨勢(shì)和規(guī)律，為疾病防控提供參考。

1 資料與方法

1.1資料來(lái)源河南省AHC月發(fā)病資料來(lái)源于河南省衛(wèi)生健康委員會(huì)(http://wsjkw.henan.gov.cn/zfxxgk/yqxx/)2013年1月至2020年12月的傳染病數(shù)據(jù)，其中2015年1月和2015年12月數(shù)據(jù)缺失，利用拉格朗日插值法對(duì)數(shù)據(jù)進(jìn)行補(bǔ)充[6]。2013年1月至2020年6月的數(shù)據(jù)用于建立模型，2020年7至12月的數(shù)據(jù)作為測(cè)試集用于驗(yàn)證模型的預(yù)測(cè)效果。

1.2建模方法SARIMA模型的一般表達(dá)式為SARIMA(p,d,q)(P,D,Q)s，其中p、P分別指非季節(jié)自回歸階數(shù)和季節(jié)自回歸階數(shù)，d、D分別指非季節(jié)差分階數(shù)和季節(jié)差分階數(shù)，q、Q分別指非季節(jié)偏自回歸階數(shù)和季節(jié)偏自回歸階數(shù)，s指季節(jié)的長(zhǎng)度[7]。采用Python3.6.4進(jìn)行建模和數(shù)據(jù)分析。具體的建模過(guò)程如下。①序列的預(yù)處理和平穩(wěn)化：將2013年1月至2020年6月的數(shù)據(jù)按月份導(dǎo)入，并以月份定義時(shí)間序列。繪制并分解時(shí)間序列，觀察時(shí)間序列的總體趨勢(shì)、季節(jié)趨勢(shì)和隨機(jī)誤差。此后對(duì)序列進(jìn)行對(duì)數(shù)轉(zhuǎn)換，并進(jìn)行非季節(jié)差分和季節(jié)差分，達(dá)到序列平穩(wěn)化的目的，同時(shí)給出d和D的值。利用迪基福勒檢驗(yàn)(augmented Dickey-fuller test,ADF)判斷該序列是否為平穩(wěn)序列。②模型的識(shí)別：繪制出平穩(wěn)序列的自相關(guān)系數(shù)(auto correlation function,ACF)圖和偏自相關(guān)系數(shù)(partial auto correlation function,PACF)圖，初步確定p、q和P、Q的值；由于P和Q一般不超過(guò)2[8]，可通過(guò)湊試法取擬合效果最佳的模型；s根據(jù)疾病的背景知識(shí)獲得。通過(guò)超參數(shù)優(yōu)化輸出平穩(wěn)序列中所有可能的p、d、q和P、D、Q的值，選擇AIC或BIC較小的模型作為備選模型。③模型的篩選與評(píng)價(jià)：首先利用Ljung-Box檢驗(yàn)進(jìn)行模型預(yù)測(cè)殘差的白噪聲檢驗(yàn)；其次通過(guò)參數(shù)值和檢驗(yàn)統(tǒng)計(jì)量，剔除參數(shù)沒(méi)有統(tǒng)計(jì)學(xué)意義的模型，篩選出最優(yōu)模型；最后對(duì)2020年7至12月的發(fā)病人數(shù)進(jìn)行擬合，計(jì)算平均絕對(duì)誤差和平均相對(duì)誤差。④模型的預(yù)測(cè)：利用最優(yōu)模型預(yù)測(cè)2021年1至12月河南省AHC發(fā)病人數(shù)。

1.3統(tǒng)計(jì)學(xué)處理采用Excel軟件建立數(shù)據(jù)庫(kù)，用Python3.6.4中“pandas”模塊和“numpy”模塊進(jìn)行數(shù)據(jù)處理，“matplotlib”模塊實(shí)現(xiàn)數(shù)據(jù)可視化，“statsmodels”模塊建立時(shí)間序列模型。檢驗(yàn)水準(zhǔn)α=0.05。

2 結(jié)果

2.1發(fā)病人數(shù)利用拉格朗日插值法估計(jì)出2015年1月的發(fā)病人數(shù)為168，2015年12月的發(fā)病人數(shù)為184。2013年1月至2020年6月河南省AHC年均發(fā)病人數(shù)為2 556，月均發(fā)病人數(shù)為213，發(fā)病人數(shù)在2017年7月最多，為380。圖1展示了AHC逐月發(fā)病人數(shù)時(shí)間和季節(jié)性趨勢(shì)分解結(jié)果。長(zhǎng)期趨勢(shì)分析結(jié)果表明，2013到2020年，發(fā)病人數(shù)總體呈先上升后平穩(wěn)的趨勢(shì)。季節(jié)分解結(jié)果顯示該病具有明顯的季節(jié)性周期分布特點(diǎn)，每年的6至8月為發(fā)病的高峰期，隨機(jī)誤差保持在一定水平范圍內(nèi)。

圖1 2013年1月至2020年6月河南省AHC發(fā)病人數(shù)的時(shí)間序列

2.2時(shí)間序列的平穩(wěn)化隨時(shí)間的增加時(shí)間序列Yt的波動(dòng)性變化較大，表現(xiàn)出明顯的異方差。對(duì)Yt進(jìn)行對(duì)數(shù)變換，消除異方差，得到對(duì)數(shù)變換后的序列l(wèi)nYt。對(duì)lnYt進(jìn)行ADF檢驗(yàn)，結(jié)果(表1)顯示，數(shù)據(jù)為非平穩(wěn)序列(P=0.400)。采用1階12步差分法消除時(shí)間序列的趨勢(shì)和季節(jié)影響后，ADF檢驗(yàn)顯示，該序列為平穩(wěn)序列(P=0.027)，符合SARIMA模型對(duì)于平穩(wěn)性的要求。

表1 lnYt和差分后lnYt的ADF檢驗(yàn)結(jié)果

2.3模型的識(shí)別根據(jù)差分結(jié)果，可確定s=12，d=D=1，模型表達(dá)式為SARIMA(p,1,q)(P,1,Q)12;觀察差分后序列的ACF圖和PACF圖(圖2)，均顯示1階拖尾，確定P=1,q=1;P、Q分別取0、1、2逐個(gè)嘗試，選擇AIC和BIC最小的模型作為備選模型，SARIMA(1,1,1)(0,1,1)12模型符合要求，納入為備選模型。此外，使用超參數(shù)優(yōu)化輸出p、q、P和Q取值不超過(guò)2的模型，除SARIMA(1,1,1)(0,1,1)12模型外，將AIC或BIC最小的模型納入為備選模型，又選出2個(gè)備選模型SARIMA(0,1,1)(0,1,1)12和SARIMA(2,1,1)(0,1,1)12。

圖2 差分后序列的ACF圖和PACF圖

2.4模型的篩選與評(píng)價(jià)對(duì)3個(gè)備選模型逐個(gè)擬合，結(jié)果見(jiàn)表2。3個(gè)模型殘差的Ljung-Box檢驗(yàn)結(jié)果顯示P均大于0.05，說(shuō)明均為白噪聲序列。僅模型SARIMA(0,1,1)(0,1,1)12的所有參數(shù)均存在統(tǒng)計(jì)學(xué)意義。綜合表2結(jié)果可知，模型SARIMA(0,1,1)(0,1,1)12的表達(dá)式為：ΔΔ12lnYt=(1-0.730L)(1-0.671L12)εt，其中L為后移算子，Δ為非季節(jié)差分算子，Δ=1-L，Δ12為12期季節(jié)差分算子，Δ12=1-L12，εt為白噪聲序列。

利用模型SARIMA(0,1,1)(0,1,1)12擬合2020年7至12月河南省AHC發(fā)病人數(shù)，結(jié)果如表3所示，平均絕對(duì)誤差為24.50，平均相對(duì)誤差為10.28%，說(shuō)明該模型擬合效果較好。該模型僅在 2020年12月的預(yù)測(cè)上存在較大的相對(duì)誤差(23.00%)，其余各月的相對(duì)誤差不超過(guò)16.00%，提示該模型具有較佳的預(yù)測(cè)性能。

表3 2020年7至12月AHC逐月發(fā)病人數(shù)預(yù)測(cè)結(jié)果

2.5模型的短期預(yù)測(cè)結(jié)果利用模型SARIMA(0,1,1)(0,1,1)12對(duì)2021年1至12月河南省AHC發(fā)病人數(shù)進(jìn)行預(yù)測(cè)，結(jié)果見(jiàn)表4及圖3。預(yù)計(jì)2021年6月份發(fā)病人數(shù)達(dá)到高峰(281)，2021年1月份發(fā)病人數(shù)最少(106)，整體的發(fā)病趨勢(shì)與2020年相比略有下降。

表4 2021年1至12月河南AHC發(fā)病人數(shù)預(yù)測(cè)結(jié)果

圖3 2021年1至12月河南AHC發(fā)病人數(shù)預(yù)測(cè)

3 討論

自2013年以來(lái)，河南省AHC的發(fā)病人數(shù)呈增長(zhǎng)趨勢(shì)，是不容忽視的公共衛(wèi)生問(wèn)題。AHC全年均可發(fā)生，每年的6至8月為發(fā)病高峰期，表現(xiàn)出明顯的季節(jié)性和周期性特征，這與先前的一些研究[9]結(jié)果相符。目前尚無(wú)用于預(yù)防AHC的疫苗或治療AHC的抗病毒藥物，但AHC通常會(huì)自行消退，不需進(jìn)一步治療[10]。河南省夏季炎熱且降雨頻繁，是AHC流行的高峰期，公共衛(wèi)生部門(mén)應(yīng)加強(qiáng)對(duì)游泳池、浴池、理發(fā)室等公共場(chǎng)所的衛(wèi)生管理與監(jiān)督，加強(qiáng)健康教育和提高居民的健康素養(yǎng)是預(yù)防和控制AHC暴發(fā)的關(guān)鍵。

SARIMA模型作為經(jīng)典的時(shí)間序列模型，在醫(yī)療資源、意外傷害、傳染病發(fā)病預(yù)測(cè)等醫(yī)療衛(wèi)生領(lǐng)域有著廣泛應(yīng)用[11-13]。SARIMA模型可以很好地獲取序列的周期性變化規(guī)律，適用于季節(jié)性或非季節(jié)性數(shù)據(jù)，具有較好的短期預(yù)測(cè)效果。本研究先通過(guò)圖示法合并湊試法確定1個(gè)SARIMA模型。由于圖示法帶有主觀性，且輸出模型的精度并非最高，故再通過(guò)超參數(shù)優(yōu)化輸出2個(gè)模型。經(jīng)過(guò)參數(shù)檢驗(yàn)、AIC、BIC和擬合效果的綜合評(píng)價(jià)，最終選出最優(yōu)模型SARIMA(0,1,1)(0,1,1)12。使用最優(yōu)模型對(duì)2020年7至12月AHC發(fā)病人數(shù)進(jìn)行擬合，平均相對(duì)誤差為10.28%，其中2020年10月及11月的數(shù)據(jù)擬合效果很好，相對(duì)誤差僅為1.09%和4.97%，反映出模型整體的預(yù)測(cè)效果較好；但2020年12月的預(yù)測(cè)值存在較大的相對(duì)誤差(23.00%)，反映出SARIMA模型不能有效提取時(shí)間序列中的隨機(jī)信息。SARIMA(0,1,1)(0,1,1)12的預(yù)測(cè)結(jié)果提示未來(lái)一段時(shí)間內(nèi)，AHC仍然是一個(gè)威脅當(dāng)?shù)厝巳航】档闹匾獋魅静　?/p>

本研究仍存在局限性。AHC病例數(shù)據(jù)的收集主要通過(guò)臨床診斷。由于居民的就診意識(shí)差、AHC自愈性強(qiáng)且預(yù)后良好，該病就診人數(shù)相對(duì)較少；此外，部分醫(yī)療機(jī)構(gòu)工作疏忽、鄉(xiāng)鎮(zhèn)醫(yī)療單位因診治能力較差導(dǎo)致AHC漏診漏報(bào)現(xiàn)象普遍存在；由于漏報(bào)的存在，預(yù)測(cè)值反映的是未來(lái)一段時(shí)間內(nèi)AHC報(bào)告例數(shù)，而非實(shí)際發(fā)病人數(shù)。此外，本研究中2015年1月和12月數(shù)據(jù)缺失。由于在對(duì)醫(yī)學(xué)時(shí)間序列數(shù)據(jù)進(jìn)行建模預(yù)測(cè)時(shí),序列的長(zhǎng)度和完整性對(duì)擬合的可靠性有影響，醫(yī)學(xué)時(shí)間序列的觀測(cè)值具有不可重復(fù)的特點(diǎn),缺失數(shù)據(jù)的插補(bǔ)或跳過(guò),會(huì)使擬合結(jié)果難以很好地反映縱向數(shù)據(jù)的規(guī)律,制約了時(shí)間序列在醫(yī)學(xué)領(lǐng)域的應(yīng)用[14]。當(dāng)缺值點(diǎn)不超過(guò)3個(gè)時(shí)，采用拉格朗日插值的效果較好[6]，故本研究中采用拉格朗日插值法進(jìn)行缺失值填補(bǔ)，但填補(bǔ)結(jié)果可能會(huì)存在一定的偏差，結(jié)果解釋需要謹(jǐn)慎。本模型在對(duì)2020年7至12月AHC發(fā)病人數(shù)進(jìn)行擬合時(shí)，各月的相對(duì)誤差大多穩(wěn)定在15%以內(nèi)，擬合效果較好，提示漏報(bào)情況在各月份中均勻分布，填補(bǔ)結(jié)果偏差較小，對(duì)模型的預(yù)測(cè)效果影響較小，仍可進(jìn)行AHC的發(fā)病趨勢(shì)預(yù)測(cè)。

綜上所述，本研究利用2013年1月至2020年6月河南省AHC發(fā)病數(shù)據(jù)(數(shù)據(jù)于刊發(fā)前更新)建立了SARIMA模型，探討了該病的發(fā)病規(guī)律并預(yù)測(cè)了2021年1至12月發(fā)病人數(shù)，該模型對(duì)河南省AHC的監(jiān)測(cè)有一定的應(yīng)用價(jià)值。由于SARIMA模型不能有效地提取時(shí)間序列中的隨機(jī)信息，難以從環(huán)境、生物、社會(huì)等因素探討疾病發(fā)生、發(fā)展和傳播特點(diǎn)，預(yù)測(cè)值可能會(huì)出現(xiàn)較大波動(dòng)。可考慮在今后的工作中優(yōu)化模型，如將SARIMA模型與支持向量機(jī)、隨機(jī)森林等模型結(jié)合，以提高預(yù)測(cè)的精度。