差分自回歸移動平均模型在南通市手足口病疫情預(yù)測中的應(yīng)用

練 維，魏 葉，韓穎穎，帥小博

1南通市疾病預(yù)防控制中心科研質(zhì)量管理科，2急性傳染病防制科，3慢性非傳染病防制科，江蘇 南通 226007；4南通市崇川區(qū)疾病預(yù)防控制中心檢驗科，江蘇 南通 226000

手足口病（hand?foot?mouth disease，HFMD）是一種病毒性的急性傳染病，1980年在我國廣州首次報道［1］，2008年被列為法定丙類傳染病，它由包括柯薩奇病毒、新型腸道病毒、埃可病毒等在內(nèi)的20 余種腸道病毒引起，其中又以腸道病毒71 型（EV71）和柯薩奇病毒A組16型（CoxA16）最為常見，一直為優(yōu)勢流行毒株，這兩者可交替或者共同傳播［2］。但近年來，由柯薩奇病毒A 組6 型（CoxA6）、柯薩奇病毒A 組9 型（CoxA9）、柯薩奇病毒A 組10 型（CoxA10）、埃可病毒11型（ECV11）等引發(fā)的疫情有上升趨勢［3-5］，值得密切關(guān)注，5歲以下尤其是3歲以下的兒童是防控的重點人群［6］。本文對2010—2019年南通市手足口病發(fā)病情況進(jìn)行分析，以2010年1月—2019年6月南通市手足口病分月報告病例數(shù)據(jù)為基礎(chǔ)，通過建立差分自回歸移動平均（autoregressive integrated moving average，ARIMA）模型［7］，預(yù)測2019 年7—12月發(fā)病率，并與當(dāng)時的實際發(fā)病率比較，驗證模型效果，為南通市手足口病的防控提供策略和技術(shù)支撐。

1 資料和方法

1.1 資料

資料來源于中國疾病預(yù)防控制信息系統(tǒng)2010年1 月1 日—2019 年12 月31 日南通市常住人口的手足口病監(jiān)測數(shù)據(jù)和人口數(shù)據(jù)。

1.2 方法

以2010年1月1日—2019年6月30日南通市手足口病疫情月報告發(fā)病率為基礎(chǔ)建立時間序列。手足口病月發(fā)病率時間序列為季節(jié)性時間序列，故采用乘積季節(jié)模型，即ARIMA（p，d，q）×（P，D，Q）S。其中d為平穩(wěn)化過程中差分的階數(shù)，p、q分別為自回歸階數(shù)和移動平均階數(shù)。P、Q分別為季節(jié)自回歸和移動平均階數(shù)，D為季節(jié)差分階數(shù)，s為季節(jié)性周期循環(huán)長度。模型建立的步驟分為：①確保時序的平穩(wěn)性。平穩(wěn)性是ARIMA模型中的一個重要假設(shè)，一般可通過時間序列圖直觀判斷。對于不平穩(wěn)時序，則需要通過數(shù)據(jù)變換和差分使序列滿足平穩(wěn)性假定，并使用ADF統(tǒng)計檢驗來驗證平穩(wěn)性假定。②模型識別和定階。根據(jù)平穩(wěn)化后序列的季節(jié)差分自相關(guān)（ACF）函數(shù)圖和偏自相關(guān)（PACF）函數(shù)圖，進(jìn)行模型的識別和定階。③模型的參數(shù)估計和檢驗。使用非線性最小二乘法估計模型的參數(shù)，Ljung?BoxQ統(tǒng)計量檢驗，P＞0.05，提示殘差是白噪聲序列，反之為非白噪聲序列。④評價模型預(yù)測效果。利用非參數(shù)檢驗法（兩個獨立樣本檢驗法），對2019年下半年實際月發(fā)病率與預(yù)測月發(fā)病率進(jìn)行比較，評價模型預(yù)測效果。

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

運用Excel2007 建立和管理數(shù)據(jù)庫，運用SPSS 25.0 軟件進(jìn)行統(tǒng)計分析和模型構(gòu)建，分類資料采用例數(shù)和構(gòu)成比進(jìn)行統(tǒng)計描述；采用非參數(shù)檢驗法對率進(jìn)行比較：利用χ2檢驗對不同型別手足口病陽性檢出率進(jìn)行比較，利用兩個獨立樣本檢驗法對手足口病實際發(fā)病率與預(yù)測發(fā)病率兩組數(shù)據(jù)進(jìn)行比較，以P＜0.05為差異有統(tǒng)計學(xué)意義。

2 結(jié)果

2.1 手足口病發(fā)病情況

2010—2019年南通市共報告手足口病90 766例，年平均發(fā)病率為124.36/10 萬，其中2018 年發(fā)病率最高，為235.29/10萬（圖1）。

圖1 2010—2019年南通市手足口病發(fā)病率Figure 1 The morbidity of HFMO in Nantong from 2010 to 2019

2.2 病原學(xué)檢測結(jié)果

2010—2019 年90 766 例中實驗室診斷病例數(shù)為2 950 例，其中EV71 陽性757 例，占25.75%；CoxA16 陽性925 例，占31.4%；其他腸道病毒陽性1 268例，占42.9%，不同的陽性型別檢出率比較，差異有統(tǒng)計學(xué)意義（χ2=206.95，P＜0.01）。不同年份間EV71、CoxA16和其他腸道病毒陽性檢出率差異均有統(tǒng)計學(xué)意義（χ2值分別為365.20、119.06、334.47，P均＜0.01，表1）。

2.3 時間分布情況

2010—2019 年，每個月份均有病例發(fā)生，經(jīng)統(tǒng)計，2014 年和2018 年發(fā)病水平較高，總體上2 月份發(fā)病數(shù)最少，6 月份發(fā)病數(shù)最多，有明顯季節(jié)性，呈雙峰特征，為夏季（5、6、7月）高峰和冬季（11、12月）次高峰，從2013年起，有隔年高發(fā)、偶數(shù)年份增強的流行特征（表2）。

表1 2010—2019年南通市手足口病病原學(xué)檢測結(jié)果Table 1 The results of etiological detection of HFMO in Nantong from 2010 to 2019

2.4 預(yù)測模型構(gòu)建

以2010年1月—2019年6月南通市手足口病分月報告病例數(shù)據(jù)為基礎(chǔ)，構(gòu)建符合季節(jié)性時間序列的ARIMA（p，d，q）×（P，D，Q）S模型，獲得2019年7—12月發(fā)病預(yù)測值，然后用2019年7—12月全市手足口病月發(fā)病率為驗證數(shù)據(jù)進(jìn)行驗證，并繪制實際值和預(yù)測值序列圖。估計預(yù)測值與實際值相對誤差來判斷模型的預(yù)測效果。

2.4.1 序列平穩(wěn)化和預(yù)測模型識別

2010年1月—2019年6月南通市手足口病月發(fā)病率時間序列圖存在以12 個月為1 個周期的季節(jié)性波動趨勢（圖2），不能滿足平穩(wěn)化的要求，需對數(shù)據(jù)進(jìn)行平穩(wěn)化處理；通過對原序列進(jìn)行一階季節(jié)差分處理后，時間序列ACF 和PACF 函數(shù)無明顯截尾和拖尾的現(xiàn)象（圖3、4），亦無線性衰減，差分后的時間序列圖接近平穩(wěn)，提示差分后序列適合時間序列模型，確定模型為ARIMA（p，d，q）×（P，D，Q）12。

2.4.2 模型的參數(shù)估計和檢驗

根據(jù)序列平穩(wěn)化處理過程，確定模型的d 和D值分別為0和1；根據(jù)手足口病的季節(jié)性調(diào)整，確定s為12；根據(jù)ACF函數(shù)圖和PACF函數(shù)圖，確定模型的p和q值均為1。模型中的P、Q值分別取0～2由低階到高階逐個進(jìn)行摸索試驗。采用Ljung?Box 方法檢驗殘差白噪聲，若為非白噪聲模型則排除。通過篩選，模型ARIMA（1，0，1）（1，1，1）12標(biāo)準(zhǔn)化的貝葉斯信息量（BIC）值為3.59，在模型中最小，R2=0.83，較大，殘差序列的自相關(guān)系數(shù)及偏自相關(guān)系數(shù)均落在95%CI 中，Ljung?BoxQ=19.77，P=0.14＞0.05，提示殘差是白噪聲序列，確認(rèn)ARMIA（1，0，1）（1，1，1）12模型最優(yōu)，可以使用。對該模型的參數(shù)進(jìn)行檢驗，差異均有統(tǒng)計學(xué)意義（P＜0.05，表3）。

表2 2010—2019年南通市手足口病各月份發(fā)病情況Table 2 Monthly incidence of HFMO in Nantong from 2010 to 2019 （例）

圖2 2010年1月—2019年6月南通市手足口病月發(fā)病率時間序列圖Figure 2 Time series plot of the monthly morbidity of HF?MO in Nantong from January 2010 to June 2019

圖3 2010年1月—2019年6月南通市手足口病月發(fā)病率季節(jié)差分自相關(guān)（ACF）函數(shù)圖Figure 3 Seasonal difference autocorrelation function graph of the monthly morbidity of HFMO in Nantong from January 2010 to June 2019

圖4 2010年1月—2019年6月南通市手足口病月發(fā)病率季節(jié)差分偏自相關(guān)（PACF）函數(shù)圖Figure 4 Seasonal difference partial autocorrelation func?tion graph of the monthly morbidity of HFMO in Nantong from January 2010 to June 2019

2.4.3 預(yù)測模型的效果

圖5 模型ARIMA（1，0，1）（1，1，1）12殘差序列的自相關(guān)與偏自相關(guān)圖Figure 5 Autocorrelation and partial autocorrelation function graphs of ARIMA（1，0，1）（1，1，1）12 model residual sequence

表3 ARIMA（1，0，1）（1，1，1）12模型參數(shù)估計Table 3 Parameter estimation of ARIMA（1，0，1）（1，1，1）12 model

通過ARIMA（1，0，1）（1，1，1）12模型對2019 年7—12 月手足口病發(fā)病率進(jìn)行預(yù)測，預(yù)測結(jié)果分別為7.08/10 萬、1.81/10 萬、3.74/10 萬、7.21/10 萬、10.71/10萬和11.29/10萬，繪制實際值和預(yù)測值序列圖（圖7）。

采用兩個獨立樣本檢驗法，比較2019 年7—12 月手足口病實際發(fā)病率（10.96/10萬、4.30/10萬、5.15/10 萬、6.25/10 萬、5.83/10 萬和3.82/10 萬）與預(yù)測值的差異，經(jīng)分析兩者間的差異無統(tǒng)計學(xué)意義（Z=0.48，P=0.63），預(yù)測數(shù)據(jù)與真實值比較接近，因此該模型有較好的預(yù)測效果。

3 討論

圖6 ARIMA（1，0，1）（1，1，1）12模型擬合圖Figure 6 Fitted figure of ARIMA（1，0，1）（1，1，1）12 model

本研究結(jié)果顯示，2010—2019 年南通市手足口病年平均發(fā)病率為124.36/10萬，低于同時期蘇州［8］、無錫［9］和全國［6］平均發(fā)病水平，整體呈現(xiàn)隔年高發(fā)、偶數(shù)年份增強的流行特點。值得關(guān)注的是相比往年，2017 年發(fā)病率有較明顯下降，但隨后2018 年出現(xiàn)最高峰，2019年又出現(xiàn)類似情況，按照本市手足口病的流行規(guī)律，2020年仍有很大概率出現(xiàn)發(fā)病高峰，這可能與易感人群的積累有關(guān)［10］。因此做到關(guān)口前移，提前做好宣傳防控并加大疫苗接種的覆蓋率很有必要。

從病原學(xué)檢測結(jié)果來看，2010—2012年引起南通市手足口病的主要病原體為EV71 或CoxA16，從2013 年起轉(zhuǎn)變?yōu)槠渌c道病毒，這種變化趨勢與國內(nèi)其他地方報道相一致［11-12］。從2010—2019年各月份發(fā)病數(shù)來看，疫情有明顯季節(jié)性，呈雙峰特征，為夏季（5、6、7月）高峰和冬季（11、12月）次高峰，亦符合我國大部分省份手足口病的發(fā)病模式［13］，究其原因，主要是5、6、7月氣溫適宜，濕度較大，有利于腸道病毒的生存和繁殖，而11、12月兒童的自身免疫力下降。

由于腸道病毒各型之間沒有交叉免疫力，且病原體類型較多，傳播途徑復(fù)雜，造成了手足口病反復(fù)感染，發(fā)病率一直居高不下。因此提前建立預(yù)警模型，進(jìn)行預(yù)警分析，進(jìn)而采取相應(yīng)措施就顯得很有必要。ARIMA模型主要用于擬合具有平穩(wěn)性（或者可以被轉(zhuǎn)換成平穩(wěn)序列）的時間序列，最早由Box和Jenkins提出，屬于時間序列模型，在手足口病、肺結(jié)核、猩紅熱等傳染病的預(yù)測中得到了廣泛的應(yīng)用，具有較高的實際價值［14］。本文以2010年1月—2019 年6 月南通市手足口病月發(fā)病數(shù)為基礎(chǔ)，通過序列平穩(wěn)化、預(yù)測模型的識別、預(yù)測模型的參數(shù)估計和檢驗，構(gòu)建了ARIMA（1，0，1）（1，1，1）12模型，最終模型的預(yù)測結(jié)果與實際擬合度較高，表明該模型對手足口病的流行趨勢預(yù)測有很好的適用性。盡管其具有擬合度較好、適用性強的特點，尤其對周期性、季節(jié)性的傳染性疾病有明顯優(yōu)勢，但也存在一定局限性，比如只能應(yīng)用于短期預(yù)測，一般不超過時間序列的10%［15］，受外部因素影響［16］等。有資料顯示，手足口病的發(fā)病與溫度、濕度等氣象因素以及人口密度呈正相關(guān)［17］，并由此衍生了空間累加非線性時空、空間貝葉斯風(fēng)險評估等一系列模型，提示在今后的實際應(yīng)用中，可將更多的影響因素納入手足口病的預(yù)警系統(tǒng)，進(jìn)一步優(yōu)化模型，從而提高模型的整體預(yù)警效果。