李 澤, 張學良
(新疆醫(yī)科大學1公共衛(wèi)生學院, 2醫(yī)學工程技術學院, 烏魯木齊 830011)
病毒性肝炎是由肝炎病毒引發(fā)的傳染性疾病,主要分為甲、乙、丙、丁、戊共五種類型。其中丙肝(hepatitis C virus,HCV)發(fā)病率逐年上升且死亡率較高,對人類健康危害較大,主要通過輸血、靜脈毒品注射、血液透析和器官移植等途徑傳播,醫(yī)源性傳播是其主要的擴散形式[1]。我國丙肝感染人數估計約3 800萬人,50%為病毒攜帶者,是全球感染人數最多的國家[2]。2005-2014年,新疆維吾爾自治區(qū)累計報告丙肝發(fā)病數83 983例,死亡數115人。本文使用新疆地區(qū)的丙肝歷史數據,結合時間序列方法中的ARIMA(Autoregressive Integrated Moving Average)乘積季節(jié)模型建立新疆地區(qū)丙肝月發(fā)病數模型,在此基礎上進行擬合和短期預測,為丙肝的防控提供一定的依據。
1.1資料來源2005-2014年新疆丙肝月發(fā)病例數來源于公共衛(wèi)生科學數據中心。
1.2研究方法
1.2.1 ARIMA季節(jié)模型 乘積季節(jié)模型考慮了時間序列的長期趨勢、循環(huán)波動、季節(jié)變化以及隨機波動之間相互影響[3],其公式簡記為ARIMA(p, d, q)×(P, D, Q)S。其中p表示自回歸階數,d表示差分階數,q表示移動平均階數,對應的參數P、D和Q分別表示季節(jié)自回歸階數、季節(jié)差分階數和季節(jié)移動平均階數。
1.2.2 ARIMA季節(jié)模型建模步驟 (1)平穩(wěn)性檢驗:常見的平穩(wěn)性檢驗有圖檢法和單位根檢驗,如ADF檢驗、DFGLS檢驗、KPSS檢驗和NP檢驗,其中ADF檢驗和KPSS檢驗運用較多[4];(2)數據變換:非平穩(wěn)序列在經過Box-Cox變換和差分處理后可轉換為平穩(wěn)序列,它是一種將倒數變換、指數變換、對數變換結合起來的變換方法[5],同時能實現方差齊性并消除異方差[6],數據變換后需重新做平穩(wěn)性檢驗;(3)純隨機性檢驗:純隨機性檢驗選用QBP或QLB統計量,當P<0.05時認為此時間序列為非白噪聲序列,說明此平穩(wěn)序列中包含值得提取的信息;(4)確定模型結構:繪制自相關圖ACF和偏自相關圖PACF,根據表1中的規(guī)則,估算模型ARIMA(p, d, q)×(P, D, Q)S中參數p、d、P和Q的范圍,從而確定候選模型;(5)估計模型參數:使用矩估計作為最大似然估計和最小二乘法迭代的初始值,并估計各個候選模型的參數;(6)模型和參數顯著性檢驗:若模型殘差通過白噪聲檢驗且滿足方差齊性,說明此最優(yōu)模型的殘差為白噪聲,否則選擇其他次優(yōu)候選模型,其次還需對模型中的參數做顯著性檢驗,如果有任何一個參數不顯著,則不再選擇此模型,而重新選擇其他候選模型再次檢驗;(7)尋找最小信息準則模型:為了選擇其中最合理的模型,還需要計算其信息準則函數值,常見有AIC、AICc、BIC、DIC、HQC,因AIC/AICc在理論上比BIC更有優(yōu)勢[7],且當樣本量足夠大時AICc會收斂于AIC[8],同時AICc更適用于時間序列模型,因此本文選用AICc作為最優(yōu)模型的評價指標;(8)模型的交叉驗證和預測:考慮到時間序列的特點,不宜采用K-fold交叉驗證,選用Hold-Out較為合適,把時序數據劃分為訓練集和驗證集,在訓練集上建立模型并估計參數,再將候選模型的預測值和驗證集進行比較從而判斷誤差,常見的擬合效果評價指標有MSE、MAPE和SMAPE。
1.3數據處理軟件使用R語言3.3.3,預測包forecast 8.0,時間序列包tseries 0.10-38,單元根檢驗包fUnitRoots 3010.78。
2.1平穩(wěn)性檢驗使用R語言繪制2005-2014年新疆地區(qū)丙肝月發(fā)病數時序圖,見圖1。對數據做ADF和KPSS平穩(wěn)性檢驗,前者P=0.231,后者P<0.01,說明該時序是非平穩(wěn)的,需要進行數據變換。
2.2數據變換為減少結果出現異方差的可能性,直接對原始數據做λ=0的Box-Cox變換,即自然對數變換。為得到季節(jié)差分和非季節(jié)差分項,對原始序列做非平穩(wěn)序列的確定性分析,圖2中可以看出明顯的季節(jié)性變化,因此需要做1階12步季節(jié)差分,即D=1,S=12。隨后對差分數據再做平穩(wěn)性檢驗,發(fā)現依然是非平穩(wěn)的,所以嘗試1階非季節(jié)差分,即d=1,檢驗后發(fā)現此時序已平穩(wěn)。
圖1 2005-2014年新疆地區(qū)丙肝月發(fā)病數時序圖
圖2 丙肝月發(fā)病數的確定性分析
2.3純隨機性檢驗采用QLB統計量進行白噪聲檢驗,差異有統計學意義(P<0.01),說明此變換后的平穩(wěn)序列不是白噪聲,序列中包含值得提取的信息。
2.4確定模型結構繪制該平穩(wěn)序列的ACF和PACF,見圖3和圖4。根據表1的判斷方法,非季節(jié)參數q可能取值0、1、2,季節(jié)參數Q可能取值為0、1,非季節(jié)參數p可能取值為0、1、2、3,季節(jié)參數P可能取值為0、1、2。因此共有3×2×4×3=72個候選模型。
圖3 平穩(wěn)序列的自相關圖ACF
2.5估計模型參數使用R語言構建了72個候選模型,每個模型的參數均會被自動估計。
2.6模型和參數顯著性檢驗對72個候選模型的殘差做統計量的白噪聲檢驗。隨后做顯著性檢驗,自由度為2005年1月-2014年6月訓練集的月數總數114減去當前候選模型的參數數量。通過計算得到72個候選模型中,有12個模型呈現顯著性。
圖4 平穩(wěn)序列的偏自相關圖PACF
2.7尋找最小信息準則模型計算上述12個模型的AICc,見表2。ARIMA(2,1,0)×(1,1,0)12即為最優(yōu)模型。圖5是該模型的殘差平方圖,可以看出沒有明顯的趨勢,并未呈現出異方差性。表3為模型的參數顯著性檢驗,顯示所有參數均顯著非零。
表2 通過顯著性檢驗的候選模型的AICc值
圖5 最優(yōu)模型的殘差平方圖
參數模型回歸系數標準誤t值Par1-0.7720.097-7.985<0.001ar2-0.2740.096-2.8630.003sar1-0.4430.096-4.610<0.001
2.8模型的交叉驗證和預測為了驗證模型ARIMA(2,1,0)×(1,1,0)12的外推能力,將2005-2014年的月時序數據劃分為兩部分,2005年1月-2014年6月的月數據為訓練集,2014年7月-2014年12月的月數據為驗證集。做Hold-Out交叉驗證,訓練集MAPE=1.44%,驗證集MAPE= 4.80%,驗證集SMAPE=2.37%,擬合與預測效果均較好,擬合情況見圖6。可以看出,模型ARIMA(2,1,0)×(1,1,0)12在驗證集上的外推能力較好。表4給出了驗證集上的誤差,平均誤差為4.67%。圖6預測部分顯示出2015年的丙肝發(fā)病總數為11 788例,略高于2014年的11 715例,預測數據見表5,發(fā)病數峰值1 154例,出現在3月。
圖6 最優(yōu)模型的擬合、驗證和預測圖
時間實際值預測值絕對誤差相對誤差2014.79271013860.082014.8967923-44-0.052014.9842857150.022014.10753764110.012014.11913973600.062014.121010956-54-0.06
表5 2015年丙肝預測月發(fā)病數
丙肝逐漸成為突出的公共衛(wèi)生問題,給社會造成了一定的經濟負擔。本研究結果顯示,2005年新疆地區(qū)丙肝病例數較少,然而在2006年之后丙肝病例數持續(xù)增加,2008年之后相對穩(wěn)定,但發(fā)病數仍緩慢上升。2004-2010年新疆地區(qū)的法定傳染病發(fā)病率中,丙肝的發(fā)病率平均為26.45/10萬,死亡率平均為0.04/10萬[9],已成為影響新疆傳染病發(fā)病率的主要原因之一。此外,新疆地區(qū)地域遼闊,各地區(qū)間的經濟發(fā)展水平、衛(wèi)生意識和習慣差距大,易導致貧窮和疾病的惡性循環(huán),因此通過現有數據尋找適用于新疆地區(qū)的丙肝預測模型,將會為丙肝的防控提供一定幫助。
時間序列主要研究事物發(fā)展和變化的規(guī)律并預測未來趨勢,ARIMA是較為常用的平穩(wěn)時間序列擬合模型。本文中針對新疆地區(qū)丙肝發(fā)病數建立了AICc最小的ARIMA乘積季節(jié)模型ARIMA(2,1,0)×(1,1,0)12用于預測新疆地區(qū)丙肝發(fā)病數。但由于ARIMA模型更加適合短期預測,在做長期預測時,最好可以更多地考慮歷史數據,從而獲得精確的預測結果。本研究結果顯示個別候選模型在驗證集上的MAPE小于最優(yōu)模型,如模型ARIMA(0,1,2)×(1,1,0)12在驗證集上的MAPE=4.71%,優(yōu)于本研究選定模型的外推能力,但是最小信息函數受到模型的極大似然函數值和模型中未知參數個數的影響,說明它充分提取了數據的信息,對數據的建模更加充分。
[1] 陳兆云, 劉繼文, 孟存仁,等. 新疆地區(qū)漢族、維吾爾族、哈薩克族丙肝患者基因型研究[J]. 新疆醫(yī)科大學學報, 2015, 38(7): 855-857.
[2] 王曉軍, 張榮珍, 胡苑笙,等. 我國病毒性肝炎流行現狀研究[J]. 疾病監(jiān)測, 2004, 19(8): 290-292.
[3] 王燕. 時間序列分析:基于R[M]. 北京: 中國人民大學出版社, 2015:158.
[4] 陳雙金. 時間序列單位根檢驗方法比較[D]. 成都: 電子科技大學, 2013.
[5] 吳劉倉, 黃麗, 戴琳. Box-Cox變換下聯合均值與方差模型的極大似然估計[J]. 統計與信息論壇, 2012, 27(5): 3-8.
[6] 崔玫意, 張玉虎, 陳秋華. Box-Cox正態(tài)分布及其在降雨極值分析中的應用[J]. 數理統計與管理, 2017, 36(1): 8-15.
[7] ANDERSON B. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach[M].2nd ed.Berlin: Springer-Verlag, 2002.
[8] BURNHAM K P, ANDERSON D. Multimodel inference: understanding AIC and BIC in model selection[J].Soc Methods Res,2004,33(2):261-304.
[9] 鄭強,王新旗,曹巖,等.2004~2010年新疆法定傳染病流行趨勢分析[J]. 疾病預防控制通報,2011,26(6):6-10.