基于SARIMA模型或SARIMA-GRNN組合模型的急性出血性結(jié)膜炎發(fā)病預測

袁 依，龐振陸，崔慶霞，王春平，李望晨

濰坊醫(yī)學院公共衛(wèi)生學院衛(wèi)生統(tǒng)計學系 山東濰坊 261053

急性出血性結(jié)膜炎(acute hemorrhagic conjunctivitis，AHC)俗稱紅眼病，主要是由腸道病毒70型(EV70)、柯薩奇A24病毒變種(CVA24V)引起的自限性疾病[1]。AHC傳染性極強，尤其是在夏秋季，農(nóng)村或?qū)W校、幼托機構常發(fā)生暴發(fā)流行[2]。張生奎等[3]在預測細菌性痢疾發(fā)病率時提出將季節(jié)差分自回歸移動平均(seasonal differential autoregressive moving average,SARIMA)法和Elman神經(jīng)網(wǎng)絡(Elman neural networks,ERNN)算法結(jié)合設計組合模型，發(fā)現(xiàn)組合模型可以更好地捕捉數(shù)據(jù)前后之間的潛在聯(lián)系，且預測精度高于單一模型。SARIMA模型通過差分運算可提取疫情數(shù)據(jù)中長期趨勢、季節(jié)變動信息，而廣義回歸神經(jīng)網(wǎng)絡(generalized regression neural network,GRNN)模型作為一種自適應學習非線性信息處理系統(tǒng)，有助于追蹤和擬合傳染病發(fā)病數(shù)的時序演化規(guī)律[4]。本研究利用我國2013年1月至2020年12月AHC發(fā)病監(jiān)測數(shù)據(jù)構建SARIMA模型及SARIMA-GRNN組合模型，評價模型預測效果，為AHC的科學防控工作提供決策依據(jù)。

1 資料與方法

1.1 資料來源2013年1月至2021年6月AHC發(fā)病資料搜集于國家衛(wèi)生健康委員會官方網(wǎng)站。其中2013年1月至2020年12月AHC逐月發(fā)病數(shù)據(jù)用于構建預測模型，2021年1 月至2021年6月AHC逐月發(fā)病數(shù)據(jù)用于模型預測效能評價。

1.2 SARIMA模型和SARIMA-GRNN組合模型的構建AHC月發(fā)病數(shù)用R3.6.1軟件中aTSA、forecast、tseries數(shù)據(jù)包進行流程處理，構建傳統(tǒng)SARIMA模型和SARIMA-GRNN組合模型。

SARIMA模型結(jié)構表達式為SARIMA(p,d,q)(P,D,Q)12，式中p、d、q分別表示季節(jié)內(nèi)自回歸AR階數(shù)、季節(jié)內(nèi)差分階數(shù)、季節(jié)內(nèi)移動平均過程MA階數(shù)。s、P、D、Q分別表示季節(jié)周期時間長度、季節(jié)間自回歸過程SAR階數(shù)、季節(jié)間差分階數(shù)、季節(jié)間移動平均過程SMA階數(shù)。首先構建AHC發(fā)病監(jiān)測數(shù)據(jù)時序分布圖，并通過單位根檢驗其平穩(wěn)性。調(diào)用nsdiffs()、ndiffs()函數(shù)輸出結(jié)果為1和0，確認只需進行一次12階季節(jié)性差分運算可以消除原始數(shù)據(jù)中的短期季節(jié)性或長期趨勢性特征。繪制經(jīng)差分運算后生成的平穩(wěn)序列自相關系數(shù)圖(ACF)和偏自相關系數(shù)圖(PACF)，觀察數(shù)據(jù)短期和季節(jié)自相關特征，根據(jù)序列特征確定d、D、s取值分別為0、1和12，p、q取值為1和0，P、Q取值0或1。據(jù)此初步納入4個備選模型，對其殘差序列進行白噪聲檢驗并估計模型參數(shù)，以AIC值最小化為準則確定最優(yōu)模型。

由于SARIMA模型的季節(jié)性差分導致了13個數(shù)據(jù)缺失，因而選取2014年2月至2020年12月全國AHC月發(fā)病數(shù)的SARIMA模型擬合值作為GRNN模型輸入層，實際月發(fā)病數(shù)作為GRNN模型輸出層，構建一維輸入、一維輸出的GRNN模型。調(diào)用mapminmax()函數(shù)對訓練數(shù)據(jù)進行建模前尺度變換預處理，使得所有數(shù)據(jù)都位于0和1之間，以消除因數(shù)據(jù)相差過大而增加預測誤差的影響。采用隨機randint( )函數(shù)在83個樣本數(shù)據(jù)中隨機選取兩個樣本數(shù)據(jù)(待估點)，剩余樣本數(shù)據(jù)訓練神經(jīng)網(wǎng)絡。通過交叉驗證法獲得使函數(shù)誤差最小的最優(yōu)光滑因子(spread)，最后借助newgrnn()函數(shù)，用歸一化后的輸入變量、輸出變量和最優(yōu)spread訓練網(wǎng)絡。

1.3 模型評價采用MATLAB 和R3.6.1軟件提取2013年1月至2020年12月全國AHC月發(fā)病數(shù)演化規(guī)律性特征，與2021年1月至2021年6月AHC月發(fā)病數(shù)比較，結(jié)合誤差指標評價模型的擬合與外推性能，誤差指標包括平均絕對誤差百分比(mean absolute error percentage，MAPE)、平均絕對誤差(mean absolute error，MAE)和均方根誤差(root mean square error，RMSE)。檢驗水準α=0.05。

2 結(jié)果

2.1 AHC發(fā)病數(shù)的分布特征2013年1月至2020年12月全國AHC月發(fā)病數(shù)的時序分布如圖1所示。從長期整體分析來看，我國AHC發(fā)病呈以12個月為周期的季節(jié)性波動特點，并隨時間呈現(xiàn)單調(diào)遞增或遞減的變化趨勢；除6至9月呈現(xiàn)季節(jié)性發(fā)病高峰，其他月份發(fā)病相對穩(wěn)定。另外，2014年9月(8 485例)、2019年7月(5 265例)分別有兩次達峰值。

圖1 AHC發(fā)病數(shù)的時序分布圖

2.2 SARIMA模型原始序列差分后的時序資料圖、ACF和PACF圖分別見圖2、3。

圖2 AHC發(fā)病數(shù)的時序資料差分效果

圖3 1步12階差分后的ACF和PACF圖

4個備選模型SARIMA(1,0,0)(0,1,0)12、SARIMA(1,0,0)(0,1,1)12、SARIMA(1,0,0)(1,1,0)12、SARIMA(1,0,0)(1,1,1)12的白噪聲檢驗與模型參數(shù)估計比較見表1，最終確定最優(yōu)模型為SARIMA(1,0,0)(0,1,1)12。

表1 SARIMA備選模型參數(shù)估計與檢驗

2.3 SARIMA-GRNN組合模型經(jīng)過反復訓練嘗試，使spread從0.01開始，每次增加0.01間隔單位一直增加到1。當spread=0.11 時，所對應的待估點預測值與實際值誤差序列的均方根(RMSE) 對應取值最小，非線性逼近能力最強，此時SARIMA-GRNN組合模型擬合預測效果最佳。

2.4 兩種模型評價模型預測結(jié)果見表2。由表2可知，SARIMA(1,0,0)(0,1,1)12預測的相對誤差范圍為5.92%～42.82%，平均相對誤差為18.21%。SARIMA-GRNN組合模型預測的相對誤差范圍為0.86%～31.41%，平均相對誤差為13.07%。

表2 2021年我國AHC月發(fā)病數(shù)預測的模型比較

對兩模型擬合效果和預測能力的評價結(jié)果見圖4和表3。由圖4、表3結(jié)果可知，SARIMA模型擬合效果優(yōu)于SARIMA-GRNN組合模型，SARIMA-GRNN組合模型擬合結(jié)果未能理想逼近發(fā)病數(shù)，卻能反映發(fā)病數(shù)的高峰波。SARIMA預測誤差較大，綜合而言，組合模型優(yōu)于傳統(tǒng)模型。

橫坐標1代表起始時間為2014年2月， 84～89代表預測時間2021年1至6月

表3 兩種模型擬合與預測效果比較

3 討論

根據(jù)2013年1月至2020年12月全國AHC月發(fā)病數(shù)據(jù)，我國AHC全年皆有發(fā)病，月發(fā)病數(shù)逐年增長且以6至9月為季節(jié)性發(fā)病的高峰階段。研究[5]表明，由于引發(fā)AHC的病毒較易發(fā)生型株變替或表面抗原變異，大約每隔3～5 a我國AHC出現(xiàn)一次發(fā)病高峰。本研究數(shù)據(jù)顯示2014年與2019年出現(xiàn)了兩次AHC發(fā)病的小高峰，與上述研究相符。

由2014年2月至2020年12月我國AHC月發(fā)病數(shù)構建SARIMA模型和SARIMA-GRNN組合模型，根據(jù)誤差評價指標MAPE、MAE和RMSE進行模型擬合、預測性能比較。結(jié)果顯示：兩種模型對于時序資料周期波動和演化趨勢特征的反映與實際情況大致吻合，提示均可用于模擬AHC發(fā)病數(shù)的演化規(guī)律，但進一步研究發(fā)現(xiàn)SARIMA-GRNN組合模型預測值的誤差評價指標均較SARIMA模型小，提示SARIMA-GRNN組合模型外推預測性能優(yōu)于傳統(tǒng)SARIMA模型，與其他流行病領域的文獻[6-8]結(jié)果一致，適用于有季節(jié)變化特征的流行性傳染病的預測，可為疾病預防控制提供定量依據(jù)。

SARIMA-GRNN組合模型作為基于傳統(tǒng)模型的改進模型，可用于我國AHC月發(fā)病數(shù)演化規(guī)律分析和預測預警。在實際工作中應注意到一方面組合模型仍然是以歷史數(shù)據(jù)隨時間演化關系為依據(jù)，若歷史資料本身的演化規(guī)律性差，必然會直接影響模型提取信息的效能；另一方面?zhèn)魅静∫资芙臃N、氣候、監(jiān)測手段等因素影響，當AHC某些時點數(shù)據(jù)受到外擾因素影響而發(fā)生突變異動時，即便模型回代內(nèi)插和短期預測效果不錯，長期預測也會隨著時間延續(xù)而精度變差。因此需要及時考慮更新和追加數(shù)據(jù)，對模型進行修正改進并擬合外推。此外，有些研究者還從時間序列分析法以外，通過對傳染病發(fā)病數(shù)影響因素篩選量化、動態(tài)數(shù)據(jù)收集，建立傳染病動力學微分方程模型[9]，以及借助計算機大數(shù)據(jù)信息平臺，結(jié)合專門機構群體研究或判斷，綜合提高疫情防控決策的科學性、時效性。