黃玉萍，傅偉杰，熊長輝，劉曉青，胡國良

人間布魯氏菌病(以下簡稱“布病”)是我國法定報(bào)告的乙類傳染病，是由布魯氏菌屬細(xì)菌引起的變態(tài)反應(yīng)性人獸共患病。經(jīng)皮膚、消化道及呼吸道傳播，臨床上以長期發(fā)熱、多汗、乏力、關(guān)節(jié)疼痛、肝脾及淋巴結(jié)腫大為特點(diǎn)，若診治不及時(shí)，不徹底，易導(dǎo)致慢性感染，可遺留骨關(guān)節(jié)不可逆的器質(zhì)性損害，使活動受限或出現(xiàn)中樞神經(jīng)系統(tǒng)后遺癥[1]，嚴(yán)重影響了人群健康及社會生產(chǎn)力，是世界上一項(xiàng)重要的公共衛(wèi)生問題[2]。2015年江西省被國家確定為南方草地開發(fā)試點(diǎn)省之一，將大大促進(jìn)江西省羊養(yǎng)殖、販運(yùn)、加工銷售、餐飲等行業(yè)的發(fā)展，羊產(chǎn)業(yè)從業(yè)人員也將迅速增加。在目前疫情形勢下，運(yùn)用時(shí)間序列模型對布病疫情進(jìn)行科學(xué)的預(yù)測，為布病防控提供科學(xué)依據(jù)顯得十分重要。

1 材料與方法

1.1材料 發(fā)病數(shù)資料來自江西省疾病預(yù)防控制中心法定傳染病報(bào)告管理系統(tǒng)，2014年1月至2017年12月布病月發(fā)病數(shù)作為擬合集建立模型，對2018年1月的發(fā)病數(shù)進(jìn)行預(yù)測。

1.2 ARIMA模型的建立

1.2.1 模型識別

1.2.1.1序列平穩(wěn)化 通過觀察原始序列或差分后序列的時(shí)間序列圖(擬合曲線在未來一段時(shí)間內(nèi)的均值和方差不發(fā)生明顯變化)、ADF檢驗(yàn)(P<0.05)獲得一個(gè)非白噪聲的平穩(wěn)序列，提取原序列趨勢效應(yīng)和季節(jié)效應(yīng)[3-4]。

1.2.1.2定階 (偏)自相關(guān)函數(shù)分別在q、p階后截尾(幾乎95%(偏)自相關(guān)系數(shù)的都落在2倍標(biāo)準(zhǔn)差范圍以內(nèi))。

1.2.2參數(shù)估計(jì) 對于平穩(wěn)序列,利用最小二乘法或最大似然法進(jìn)行參數(shù)估計(jì)，參數(shù)估計(jì)值的顯著性檢驗(yàn)可以得出模型中的一些項(xiàng)是否需要，一般要求每個(gè)參數(shù)都有統(tǒng)計(jì)學(xué)意義。如果擬合的多個(gè)模型均通過參數(shù)顯著性檢驗(yàn)，引入赤池信息準(zhǔn)則( Akaikes information criterion，AIC)概念，進(jìn)行模型優(yōu)化，以該數(shù)值小者為優(yōu)[5-6]。

1.2.3模型檢驗(yàn) 模型的殘差序列自相關(guān)函數(shù)和偏相關(guān)函數(shù)均在置信區(qū)間內(nèi)，且P>0.05，該序列通過殘差白噪聲檢驗(yàn)，模型可用來預(yù)測。

1.3模型預(yù)測 模型預(yù)測采用均方根誤差(Rood Mean Squared Error， RMSE)、平均絕對誤差(Mean Absolute Error，MAE)、希爾不等系數(shù)(Theil Inequality Coefficient，Theil IC)、偏方差(Bias Proportion，BP)、方差率(Variance Proportion，VP)、協(xié)變率(Covariance Proportion，CP)等指標(biāo)衡量，前5個(gè)指標(biāo)越小、最后一個(gè)指標(biāo)越大說明預(yù)測效果越好。模型的建立和預(yù)測通過Eviews10軟件實(shí)現(xiàn)。點(diǎn)預(yù)測采用誤差進(jìn)行預(yù)測準(zhǔn)確性的比較。

2 結(jié) 果

2.1江西省2014-2017年布病月發(fā)病數(shù)序列圖 2014-2017年間江西省布病發(fā)病總數(shù)218例，單月最多12例，最少0例，數(shù)據(jù)整體不算平穩(wěn)(見圖1)。(偏)自相關(guān)圖和ADF檢驗(yàn)(P<0.05)均顯示序列平穩(wěn)(見圖2)。用Eviews10對數(shù)據(jù)進(jìn)行擬合,得到不做差分時(shí)的最優(yōu)模型AR(1)，R2=0.23，AIC=5.138，RMSE=2.934，MAE=2.429，Theil IC=0.279，BP=0.002，VP=0.339，CP=0.659。

圖1 江西省2014-2017年布病月發(fā)病數(shù)序列圖Fig.1 Monthly incidence of brucellosis in Jiangxi Province,2014 to 2017

圖2 江西省2014-2017年布病月發(fā)病數(shù)序列的(偏)自相關(guān)圖Fig.2 ACF and PACF diagram of the number of brucellosis cases per month in Jiangxi Province,2014 to 2017

2.2一階差分分析 用Eviews10對一階差分后的數(shù)據(jù)進(jìn)行擬合，得到最優(yōu)模型ARIMA(1,1,3)，R2=0.297，AIC=5.204，RMSE=2.933，MAE=2.376，Theil IC=0.282，BP=0.053，VP=0.145，CP=0.801。

2.3二階差分分析 用Eviews10對二階差分后的數(shù)據(jù)進(jìn)行擬合，得到最優(yōu)模型ARIMA(3,2,0)，R2=0.453，AIC=5.776，RMSE=4.065，MAE=3.089，Theil IC=0.326，BP=0.0003，VP=0.044，CP=0.955。

2.4三階差分分析 用Eviews10對三階差分后的數(shù)據(jù)進(jìn)行擬合，得到最優(yōu)模型ARIMA(3,3,0)，R2=0.690，AIC=6.335，RMSE=5.299，MAE=4.247，Theil IC=0.384，BP=0.0003，VP=0.253，CP=0.747。

2.5模型的擬合度和預(yù)測誤差分析 從擬合圖(見圖3)來看，擬合曲線最吻合的為ARIMA(3,3,0)，其次為ARIMA(3,2,0)；AR(1)和ARIMA(1,1,3)擬合效果較差。就各模型擬合的指標(biāo)(見表1)來看，AR(1)的AIC值(5.138)和希爾不等指數(shù)(0.279)最小；ARIMA(1,1,3)的均方根誤差(2.933)和平均絕對誤差(2.376)最??；ARIMA(3,2,0)的偏差率(0.0003)和方差率(0.044)最小，協(xié)變率(0.955)最大；ARIMA(3,3,0)的R2值(0.69)最大。

2018年1月江西省布病報(bào)告數(shù)為5例，4個(gè)最優(yōu)模型中，AR(1)的預(yù)測誤差最小為0.21，其次是ARIMA(3,3,0)。拋開預(yù)測誤差較大的ARIMA(1,1,3)和ARIMA(3,2,0)，比較AR(1)和ARIMA(3,3,0)，可以發(fā)現(xiàn)ARIMA(3,3,0)有著更吻合的趨勢圖和最大的R2值，而AR(1)有著更小的AIC值和希爾不等指數(shù)。

圖3 四種模型擬合圖的比較Fig.3 Comparison of fitting graphs of four models

表1 四種模型擬合參數(shù)的比較
Tab.1 Comparison of fitting indicators of four models

模型模型擬合指標(biāo)預(yù)測值誤差R2AICRMSEMAETheil ICBPVPCPAR(1)0.2305.1382.9342.4290.2790.0020.3390.6595.210.21ARIMA(1,1,3)0.2975.2042.9332.3760.2820.0530.1450.8016.301.3ARIMA(3,2,0)0.4535.7764.0653.0890.3260.00030.0440.9556.371.37ARIMA(3,3,0)0.6906.3355.2994.2470.3840.00030.2530.7474.10-0.9

3 討 論

江西省2014-2017年布病報(bào)告總數(shù)為218例，病例多集中于5-9月，與北方地區(qū)不一致[7]。發(fā)病高峰不穩(wěn)定。原因可能是布病近年來才在江西省內(nèi)出現(xiàn)，相關(guān)醫(yī)務(wù)工作者及患者對其較陌生，對于發(fā)現(xiàn)疾病的靈敏度不夠，從而出現(xiàn)了診斷和報(bào)告的滯后現(xiàn)象。往年1月份江西省布病報(bào)告數(shù)一般為0～2例，2018年1月報(bào)告數(shù)為5例，明顯超過以往水平。報(bào)告單位都是省內(nèi)市級及以上醫(yī)療單位，2017年底江西省內(nèi)布病病例分布覆蓋所有地級市[8]，這說明各市級醫(yī)療單位對于布病有了一定的警惕性，但同時(shí)要加強(qiáng)縣(區(qū))級醫(yī)療機(jī)構(gòu)的人員培訓(xùn)，以實(shí)現(xiàn)布病的早發(fā)現(xiàn)、早治療，早處理、早報(bào)告。值得關(guān)注的是，流調(diào)報(bào)告顯示該5例病例中，3例有明顯的羊接觸史；1例為1歲兒童，食用新鮮蒸煮過的羊奶；1例為家務(wù)待業(yè)者，無明確傳染源。結(jié)合病例中非職業(yè)人群的占比由2015年的30%增長至2017年的48%，表明布病的傳播形式多樣，疾病的溯源和防控依然存在著較大難度。

ARIMA模型是根據(jù)數(shù)據(jù)序列偏、自相關(guān)函數(shù)建立起線性的數(shù)據(jù)間相互依賴的定量模型，整合了趨勢因素、周期因素和隨機(jī)誤差等因素，既吸取了回歸分析的優(yōu)點(diǎn)又發(fā)揮了移動平均的長處，且只需發(fā)病資料就可進(jìn)行預(yù)測，是短期預(yù)測最為成熟的時(shí)間序列預(yù)測方法之一[9-10]。

本研究探索在不同差分次數(shù)條件下創(chuàng)建預(yù)測模型，得出了幾個(gè)最優(yōu)模型。其中，AR(1)和ARIMA(3,3,0)的預(yù)測值最接近于實(shí)際值，而ARIMA(3,2,0)和ARIMA(1,1,3)的預(yù)測誤差較大，說明赤池信息準(zhǔn)則和R2值是確定最優(yōu)模型的兩個(gè)較好指標(biāo)。另一方面，ARIMA(3,2,0)的R2、協(xié)方率和擬合圖的表現(xiàn)都好于ARIMA(1,1,3)，預(yù)測值的誤差卻較大，猜測可能是由于預(yù)測誤差本身存在的影響。同時(shí)，一般認(rèn)為ARIMA模型法至少需要50個(gè)以上數(shù)據(jù)建模[11]，本研究只有48個(gè)數(shù)據(jù)且含有3個(gè)零值的小樣本。

模型的選擇很大程度上是研究者根據(jù)ACF圖和PACF圖再結(jié)合自身經(jīng)驗(yàn)綜合確定的相對最優(yōu)模型，所以研究者給出的并不是最優(yōu)模型而是基于自身經(jīng)驗(yàn)做出的最優(yōu)選擇[12]。專家指出p、q的值一般不超過5，然而有研究跳出這個(gè)限制嘗試建立了幾個(gè)可能的模型，發(fā)現(xiàn)ARIMA(7,1,0)(1,0,1)12有著更小的絕對誤差和相對誤差(平均值<5%)，并且殘差檢驗(yàn)大于0.8提示擬合效果更好[13]。盡管ARIMA模型在很大程度上提取了時(shí)間序列信息，但原理上依然是一個(gè)線性模型。而布病的發(fā)生受不同宿主動物與傳播途徑等因素的影響，因此建立的預(yù)測模型并不唯一，在今后的研究中考慮在條件允許的情況下盡可能的全面收集布病發(fā)病的影響因素, 將相關(guān)因素納入模型中提高模型精度，以達(dá)到最佳的預(yù)測效果。

利益沖突：無