有序聚類分析及周期圖法在猩紅熱流行周期中的應(yīng)用研究*

段 禹黃曉磊王玉杰張俊青蘇 虹潘海峰王 靜△

·論著·

有序聚類分析及周期圖法在猩紅熱流行周期中的應(yīng)用研究*

段 禹1黃曉磊1王玉杰1張俊青2蘇 虹1潘海峰1王 靜1△

目的探討合肥市猩紅熱的流行周期，為猩紅熱發(fā)病預(yù)測和早期預(yù)警提供理論基礎(chǔ)。方法合肥市疾病預(yù)防控制中心提供1985－2003年猩紅熱病例資料的監(jiān)測數(shù)據(jù)及2004－2008年網(wǎng)絡(luò)直報監(jiān)測數(shù)據(jù)。使用有序聚類分析對猩紅熱發(fā)病階段劃分類別，使用周期圖法提取潛在周期并建立相應(yīng)的周期函數(shù)擬合發(fā)病率資料。結(jié)果1985－2008年合肥市共有1996名猩紅熱病例，年平均發(fā)病率為1．9620/10萬。24年猩紅熱發(fā)病率波動總體可按有序聚類分為低－高－低－高4個階段，分別為1985－1988年，1989－1997年，1998－2003年，2004－2008年。其中前3階段總和為19年，與周期圖法檢測出猩紅熱發(fā)病率序列存在的第一隱含周期T1＝19相同，此外序列還包含第二隱含周期T2＝5。結(jié)論有序聚類分析和周期圖法可以運(yùn)用于猩紅熱流行周期的識別和提取。

猩紅熱 流行周期 周期圖法 有序樣品聚類分析

猩紅熱是一種急性呼吸道傳染病，由A組β型鏈球菌引起，為我國法定報告的乙類傳染?。?］。猩紅熱感染者多為兒童、青少年［2］，引起的主要癥狀為發(fā)熱、咽峽炎、全身彌漫性鮮紅色皮疹和皮疹消退后明顯脫屑［3］。目前，該病尚無特異性預(yù)防疫苗，同時抗生素濫用引起的鏈球菌耐藥也為該病防治帶來挑戰(zhàn)［4］。了解猩紅熱的流行規(guī)律，預(yù)測其發(fā)病趨勢對于控制該疾病的傳播有著重要的作用。許多研究中已經(jīng)應(yīng)用統(tǒng)計學(xué)模型對傳染病流行周期進(jìn)行研究［5－6］，如：Lima等使用小波分析技術(shù)對于百日咳在智利的流行周期進(jìn)行檢測。本次研究旨在應(yīng)用有序樣品聚類及周期圖法兩種方式對合肥市1985至2008年猩紅熱流行規(guī)律進(jìn)行分析，探討其變化趨勢及流行周期，并為猩紅熱發(fā)病預(yù)測，早期預(yù)警提供方法。

資料與方法

1．猩紅熱發(fā)病資料：合肥市疾病預(yù)防控制中心提供1985－2003年猩紅熱病例資料的監(jiān)測數(shù)據(jù)及2004－2008年網(wǎng)絡(luò)直報監(jiān)測數(shù)據(jù)。猩紅熱病例的診斷按照國家衛(wèi)生部頒布標(biāo)準(zhǔn)［7］，均為確診病例，同時滿足以下3條標(biāo)準(zhǔn)：（1）具有猩紅熱臨床癥狀表現(xiàn)；（2）咽拭子或病灶分泌物血清學(xué)分群鑒定為A組β型鏈球菌；（3）猩紅熱相關(guān)的其他實(shí)驗(yàn)室檢查出現(xiàn)陽性結(jié)果或具有可疑接觸A組β型鏈球菌的流行病學(xué)史。合肥市各年的平均人口數(shù)由國家統(tǒng)計局獲取，用于計算各年度的發(fā)病率（/10萬）。

2．有序樣品聚類：有序聚類算法是針對有序樣本的一種統(tǒng)計分類方法。它將資料按照原有次序分為若干類別，屬于特殊的條件系統(tǒng)聚類。其計算方法如下［8］：

（1）定義類的直徑

長度為n的原始序列y中包含樣本｛X（1），X（2），…，X（n－1），X（n）｝。將其按原有次序分為若干類別，設(shè)其中某一類G包含的樣本有｛X（i），X（i+1），…，X（j）｝（j＞i），記為G＝｛i，i+1，…，j｝，則類別G的均值向量為：

定義類別G直徑為：

直徑D（i，j）表示類別G內(nèi)共j－i+1個樣本的總差異，其指標(biāo)是離均差平方和。

（2）定義分類損失函數(shù)

用b（n，k）表示將n個有序樣品分為k類的某一種分法，其中分割點(diǎn)分別為i1（i1＝1），i2，i3，…，ik，定義上述分類法的損失函數(shù)為：

損失函數(shù)L［b（n，k）］即為k個類別離均差平方和的總和，其越小表示該分類方式越優(yōu)。當(dāng)一種b（n，k）使分類損失函數(shù)L最小，將該分類法記為P（n，k）。

在將n個樣品聚為k（1＜k＜n）類的過程中，首先通過上述公式（3）計算得到k＝2時最小損失函數(shù)L［P（n，2）］和此時的分割點(diǎn)，在此基礎(chǔ)上增加1個新的最優(yōu)分割點(diǎn)并計算出L［P（n，3）］，迭代該計算過程可以分別得到將n個樣品聚為k（1＜k＜n）類的最優(yōu)分法，以L［P（n，k）］為縱坐標(biāo)，k為橫坐標(biāo)，畫出損失函數(shù)變化趨勢圖，選擇出最為合適的分類數(shù)。

3．周期圖法：周期圖法是一種使用試驗(yàn)周期配合實(shí)際序列，從而找出隱含周期的方法。若實(shí)際序列中確實(shí)存在頻率為ωi的隱含周期，則周期圖IN（ωi）在周期頻率ωi處有較大的峰值，故我們可以借此判斷序列中的潛在周期［9－10］。其計算方法如下：

（1）消除原始序列y的線性趨勢

若長度為n原始序列y存在趨勢或序列均值不為0，則需首先去除趨勢成分H（t），調(diào)整后的序列為y1＝y(tǒng)－H（t）。

（2）計算傅里葉系數(shù)并檢測隱含周期

其中t為序列中各項(xiàng)的期數(shù)，τ為試驗(yàn)周期，其取值范圍為［1，n－1］的正整數(shù)，K為滿足K×τ＜n的最大正整數(shù)，當(dāng)τ使得達(dá)到最大時，t即為檢測出的隱含周期Ti，需檢驗(yàn)是否為周期震動的極大值（零假設(shè)：無周期震動）：

本次研究中取α＝0．05，其對應(yīng)的J界值為2．996。若Ji小于等于界值，零假設(shè)成立，認(rèn)為無周期震動，計算結(jié)束；若拒絕零假設(shè)，則認(rèn)為T為相應(yīng)的隱含周期，此時初相位：

振幅：

則周期為T的周期函數(shù)為：

（3）擬合周期函數(shù)

將序列y1代入步驟（2），計算出第一隱含周期及其周期函數(shù)C1（t）。隨后，令y2＝y(tǒng)1－C1（t），并重復(fù)步驟（2），進(jìn)一步提取第二隱含周期。同理可得yi＝y(tǒng)i－1－Ci－1（t），迭代該過程，直到J檢驗(yàn)結(jié)果提示序列中不再存在周期震動。預(yù)測值Y（t）的95%置信區(qū)間采用bootstrap法進(jìn)行估計［11］。對于原始樣本y重復(fù)1000次樣本含量為24的有放回抽樣，從而得到1000個bootstrap樣本，使用它們重新擬合上述線性函數(shù)和正弦函數(shù)模型，并利用模型計算相應(yīng)的預(yù)測值。綜合1000個bootstrap樣本中的各期預(yù)測值的2．5分位數(shù)和97．5分位數(shù)作為Y（t）的95%置信區(qū)間的上下限。

4．統(tǒng)計學(xué)分析

使用樣品有序聚類對于猩紅熱流行階段進(jìn)行分割；使用周期圖法建立相應(yīng)的周期函數(shù)擬合發(fā)病率資料，計算潛在周期；使用單位根檢驗(yàn)（augmented dickey-fuller test，ADF）對猩紅熱發(fā)病序列進(jìn)行平穩(wěn)性檢驗(yàn)。猩紅熱發(fā)病資料的計算處理，描述性統(tǒng)計分析，樣品有序聚類與周期圖法均使用MATLAB（version 7．0）編寫程序?qū)崿F(xiàn)。本次研究中檢驗(yàn)水準(zhǔn)α＝0．05。

結(jié) 果

1．猩紅熱發(fā)病率

合肥市24年共出現(xiàn)1996名猩紅熱病例，年平均發(fā)病率為1．9620/10萬。2008年發(fā)病率最高，達(dá)到3．8317/10萬，1987年發(fā)病率最低，達(dá)到0．5316/10萬。發(fā)病率峰值分別出現(xiàn)在1991年、1995年、2005年、2008年，并依次增高。年發(fā)病率變化相對穩(wěn)定，但自2003年以后總體呈現(xiàn)上升趨勢（圖1）。

2．1985－2008年基于猩紅熱發(fā)病率的有序聚類

將24年的猩紅熱發(fā)病率進(jìn)行有序聚類，分別計算出將其聚為2～23類的最小損失函數(shù)，損失函數(shù)趨勢圖見圖2。從圖2可以看出，損失函數(shù)在4類時出現(xiàn)折點(diǎn)，隨后下降趨勢逐漸平穩(wěn)，因此首先考慮將其聚為4類。24年聚成4類，損失函數(shù)值最小為5．7738，此時分割點(diǎn)分別為i1＝1，i2＝5，i3＝14，i4＝20。按照上述最優(yōu)分法，24年可以劃分為1985－1988年，1989－1997年，1998－2003年，2004－2008年四個階段，其他具體分類信息可見表1。

圖1 1985－2008年合肥市猩紅熱年發(fā)病率序列圖（/105）

圖2 1985－2008年合肥市猩紅熱發(fā)病率聚類最小損失函數(shù)隨分組數(shù)量變化趨勢圖

表1 合肥市1985－2008年基于猩紅熱發(fā)病率的有序聚類分組

3．周期圖法

從圖1可以看出序列整體具有上升趨勢，同時ADF檢驗(yàn)結(jié)果顯示t＝－0．2976＜－1．9507，P＝0．532，接受零假設(shè)，認(rèn)為序列不平穩(wěn)。因此使用detrend函數(shù)，去除原始序列中的線性趨勢H（t），獲得去除趨勢后的發(fā)病率序列y1并再次進(jìn)行檢驗(yàn)（t＝－2．7628＞－1．9507，P＜0．01），顯示y1序列平穩(wěn)。利用y1＝y(tǒng)－H（t）求出線性趨勢函數(shù)H（t）＝1．2306 +0．0585t。

隨后，令y2＝y(tǒng)1－C1（t），并對y2再進(jìn)行周期圖分析，得出3．4100＞2．996，P＜0．05，第二隱含周期為5，其對應(yīng)的周期函數(shù)為：

再次從y2中剔除周期函數(shù)C2（t）后得到y(tǒng)3。但對序列的周期檢驗(yàn)結(jié)果顯示J＝1．8496，P＞0．05，認(rèn)為序列y3中無周期震動，故終止序列擬合。此時，將趨勢項(xiàng)與兩個周期函數(shù)線性相加得到序列的擬合函數(shù)Y（t），Y（t）與原始發(fā)病率序列的比較情況可見圖3。

圖3 1985－2008年合肥市猩紅熱發(fā)病率序列與周期圖法擬合函數(shù)圖

4．1985－2008年合肥市猩紅熱流行特點(diǎn)

周期圖法的結(jié)果顯示合肥市猩紅熱發(fā)病率波動的第一隱含周期為19年，這與有序聚類對其階段的劃分具有相同之處：第一階段1985－1988年為周期函數(shù)C1（t）從低點(diǎn)開始上升至0的過程，隨后第二階段1989－1997年則是周期函數(shù)C1（t）位于x軸上方的階段，第三階段1998－2003年則是周期函數(shù)C1（t）處于極小值前后的年份。隨后，周期函數(shù)C1（t）重新進(jìn)入上升階段，而此時則是有序聚類結(jié)果提示的第四階段2004－2008年。前3個階段總和為19年，與第一隱含周期相同。猩紅熱4個階段的平均發(fā)病率呈現(xiàn)低－高－低－高的過程，與第一隱含周期的周期函數(shù)所展示的趨勢相同。

討 論

本次研究使用有序聚類分析及周期圖法探討了合肥市1985－2008年猩紅熱流行周期。猩紅熱年發(fā)病率整體波動于0．5316～3．8317/10萬，趨勢函數(shù)H（t）提示在此期間發(fā)病率總體具有緩慢上升的趨勢。此外發(fā)病率波動的第一隱含周期為19年，但尚未有其他研究支持這一發(fā)現(xiàn)；第二隱含周期為5年，這一周期與國內(nèi)多數(shù)猩紅熱發(fā)病規(guī)律研究的結(jié)論相同或相似［12－14］。值得注意地是自2004年，雖然發(fā)病率曲線在趨勢上與擬合結(jié)果相近，但實(shí)際數(shù)值卻明顯較高。對于這一現(xiàn)象，我們提出兩種可能的原因：一方面可能是由于猩紅熱疾病自身特點(diǎn)（如病原體的基因型）與環(huán)境因素變化而導(dǎo)致［1，15］；另一方面可能與疾病監(jiān)測力度有關(guān)，由于2003年SARS的爆發(fā)，傳染病疫情引起了各部門的廣泛關(guān)注，2004年新的中華人民共和國傳染病防治法修訂通過并開始執(zhí)行，同時也開始進(jìn)行猩紅熱疫情的網(wǎng)絡(luò)直報［16］，在這種情況下，對于猩紅熱監(jiān)測重視程度也會相應(yīng)提高，可能會增加對于該疾病的正確識別率與報告率，間接地增高了監(jiān)測數(shù)據(jù)中的發(fā)病率。

我們首先對于猩紅熱發(fā)病率在時間上進(jìn)行有序聚類分析，將24年的發(fā)病變化在時間上分割為4個階段，進(jìn)而獲取到各階段發(fā)病率的特點(diǎn)，為提取流行周期提供幫助。有序聚類分析作為一種非監(jiān)督學(xué)習(xí)的分類方法，對于各種有序資料（包括時間序列）均可按照其數(shù)據(jù)結(jié)構(gòu)特征進(jìn)行最優(yōu)劃分［17］。但同時這種算法目的主要在于將資料合理地分為若干類別以待進(jìn)一步的分析擬合，故其應(yīng)用的缺陷是不能直接得出關(guān)于序列特征的結(jié)論。因此，我們在此基礎(chǔ)上，進(jìn)一步使用周期圖法對于猩紅熱流行周期進(jìn)行推斷，并將聚類分析的結(jié)果與隱含周期相結(jié)合。目前對于傳染病時間序列研究中常用的是較為形象與直觀的時域分析，如ARIMA模型，但主要局限性在于容易出現(xiàn)過度差分與過度擬合［18－19］。而周期圖法從頻域分析的角度通過傅里葉轉(zhuǎn)換的方式提取序列中的隱含周期，并通過若干周期函數(shù)的線性相加擬合原始序列，適用于平穩(wěn)序列的周期提取。應(yīng)用于非平穩(wěn)序列時，可通過剔除線性趨勢，序列轉(zhuǎn)換等方式首先將原始序列轉(zhuǎn)化為平穩(wěn)序列。

在有序聚類法劃分出的發(fā)病率變化的第4個階段中，按照擬合周期函數(shù)的變化趨勢，此時處于上升階段，這一趨勢與實(shí)際觀測值在趨勢上基本相同，且按照其周期趨勢仍有3～5年的波動上升期。根據(jù)其他研究結(jié)果顯示，2011年前后猩紅熱在合肥［20］乃至全國其他地區(qū)［3，21］均有不同程度的發(fā)病升高。因此，我們認(rèn)為有序樣品聚類與周期圖法可以較好地用于傳染病流行階段的劃分及流行周期的檢測。最后，需要注意本次研究中第一隱含周期為19年的結(jié)果，也有可能受到諸多因素的影響，如選取序列的長短、地域的不同等，故有待于將來對于更完整的發(fā)病率序列進(jìn)行分析，做進(jìn)一步驗(yàn)證。

［1］彭曉旻，楊鵬，吳雙勝，等．北京地區(qū)2011－2014年致兒童猩紅熱A組鏈球菌emm基因型別變化特征分析．中華流行病學(xué)雜志，2015，36（12）：1397-1400．

［2］Wong SS，Yuen KY．Streptococcus pyogenes and re-emergence of scarlet fever as a public health problem．Emerg M icrobes Infect，2012，1（7）：e2．

［3］You YH，Song YY，Yan XM，et al．Molecular epidem iological characteristics of Streptococcus pyogenes strains involved in an outbreak of scarlet fever in China，2011．Biomed Environ Sci，2013，26（11）：877-885．

［4］Chen YY，Huang CT，Yao SM，et al．Molecular epidemiology of group A streptococcus causing scarlet fever in northern Taiwan，2001-2002．Diagn M icrobiol Infect Dis，2007，58（3）：289-295．

［5］Lima M，Estay SA，F(xiàn)uentes R，et al．Whooping cough dynam ics in Chile（1932-2010）：disease temporal fluctuations across a north-south gradient．BMC Infect Dis，2015，15（1）：590．

［6］申銅倩，劉文東，胡建利，等．x-11-ARIMA過程在痢疾疫情預(yù)測中的應(yīng)用研究．中國衛(wèi)生統(tǒng)計，2016，31（3）：395-398．

［7］中華人民共和國衛(wèi)生部．WS282-2008猩紅熱診斷標(biāo)準(zhǔn)．北京：人民衛(wèi)生出版社，2008．

［8］Peng Z，Bao C，Zhao Y，et al．Weighted Markov chains for forecasting and analysis in Incidence of infectious diseases in jiangsu Province，China．JBiomed Res，2010，24（3）：207-214．

［9］（美）Brockwell PJ，Davis RA．著．時間序列的理論與方法（第2版）．田錚譯．高等教育出版社，2001：257-266．

［10］喬小妮，李豐森，牛剛，等．基于周期圖法的醫(yī)院門診流量管理研究．中國數(shù)字醫(yī)學(xué)，2015，10（6）：77-79．

［11］陳峰，陸守曾，楊珉．Bootstrap估計及其應(yīng)用．中國衛(wèi)生統(tǒng)計，1997，14（5）：5-7．

［12］徐斌，黃夏萍，覃曲波．南寧市1965-2004年猩紅熱流行特征分析．實(shí)用預(yù)防醫(yī)學(xué)，2006，13（5）：1208-1210．

［13］馬昭君，營亮．2004-2013年連云港市猩紅熱流行的特征．職業(yè)與健康，2015，31（3）：348-350．

［14］周雨．1997-2006年沈陽市和平區(qū)猩紅熱資料分析．預(yù)防醫(yī)學(xué)論壇，2009，15（4）：358-359．

［15］Liang Y，Liu X，Chang H，eta1．Epidemiological andmolecular characteristics of clinical isolates of Streptococcus pyogenes collected between 2005 and 2008 from Chinese children．JMed M icrobiol，2012，61（Pt7）：975-983．

［16］李雷雷，蔣希宏，隋霞，等．中國2005-2011年猩紅熱疫情流行病學(xué)分析．中國公共衛(wèi)生，2012，28（6）：826-827．

［17］楊毅，趙國浩，秦愛民．面板數(shù)據(jù)的有序聚類分析及其應(yīng)用-以全球氣候變化聚類分析為例．統(tǒng)計與信息論壇，2012，27（7）：13-18．

［18］Zhang T，Yang M，Xiao X，etal．Spectral analysis based on fast Fourier transformation（FFT）of surveillance data：the case of scarlet fever in China．Epidemiol Infect，2014，142（3）：520-529．

［19］Chen B，Sumi A，Toyoda S，etal．Time seriesanalysisof reported cases of hand，foot，and mouth disease from 2010 to 2013 in Wuhan，China．BMC Infect Dis，2015，15（1）：495．

［20］秦薇子，張笑嫣，李萍．2007-2011年度安徽省某三級甲等醫(yī)院法定傳染病疾病譜分析．中華疾病控制雜志，2013，17（3）：251-253．

［21］Lau EH，Nishiura H，Cow ling BJ，et al．Scarlet fever outbreak，Hong Kong，2011．Emerg Infect Dis，2012，18（10）：1700-1702．

（責(zé)任編輯：劉 壯）

Application of Sequential Cluster Analysis and Periodogram M ethod in Epidem ic Trend Analysis of Scarlet Fever

Duan Yu，Huang Xiaolei，Wang Yujie，etal
（Departmentof Epidemiology and Biostatistics，School of Public Health，AnhuiMedical University（230032），Hefei）

ObjectiveWe am id to analyze epidem ic trend of scarlet fever in Hefei city and provide predictivemethods for early warning of scarlet fever．MethodsSurveillance data of scarlet fever from 1985 to 2008 were collected from centers for disease control and prevention of Hefei city．Sequential cluster analysiswas used to divide these years into several periods．Periodogram method was used to extract potential cycle and fit the time series of scarlet fever．ResultsThere were altogether 1996 cases of scarlet fever in Hefei city from 1985 to 2008．The average incidence of scarlet fever was 1．9620 per 105．During these years，four clusters were classified by sequential cluster analysis which were 1985－1988，1989－1997，1998－2003，2004－2008，respectively．The first three clusterswere totally 19 yearswhich was equal to the first potential cycle T1 of scarlet fever．In addition，the second potential cycle T2 was equal to 5 in incidence series．ConclusionSequential cluster analysis and periodogram method could be used to extract epidem ic cycles of scarlet fever incidence．

Scarlet fever；Epidem ic cycle；Periodgram method；Sequential cluster analysis

安徽省自然科學(xué)基金（編號：1408085MH159）

1．安徽醫(yī)科大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計學(xué)系（230032）

2．安徽省合肥市疾病預(yù)防控制中心

△通信作者：王靜，E-mail：jwang2006＠126．com