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

·論著·

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

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

目的探討合肥市猩紅熱的流行周期，為猩紅熱發(fā)病預測和早期預警提供理論基礎。方法合肥市疾病預防控制中心提供1985－2003年猩紅熱病例資料的監(jiān)測數據及2004－2008年網絡直報監(jiān)測數據。使用有序聚類分析對猩紅熱發(fā)病階段劃分類別，使用周期圖法提取潛在周期并建立相應的周期函數擬合發(fā)病率資料。結果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。結論有序聚類分析和周期圖法可以運用于猩紅熱流行周期的識別和提取。

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

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

資料與方法

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

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

（1）定義類的直徑

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

定義類別G直徑為：

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

（2）定義分類損失函數

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

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

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

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

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

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

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

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

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

振幅：

則周期為T的周期函數為：

（3）擬合周期函數

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

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

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

結 果

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

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

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

將24年的猩紅熱發(fā)病率進行有序聚類，分別計算出將其聚為2～23類的最小損失函數，損失函數趨勢圖見圖2。從圖2可以看出，損失函數在4類時出現折點，隨后下降趨勢逐漸平穩(wěn)，因此首先考慮將其聚為4類。24年聚成4類，損失函數值最小為5．7738，此時分割點分別為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ā)病率聚類最小損失函數隨分組數量變化趨勢圖

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

3．周期圖法

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

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

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

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

4．1985－2008年合肥市猩紅熱流行特點

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

討 論

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

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

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

［1］彭曉旻，楊鵬，吳雙勝，等．北京地區(qū)2011－2014年致兒童猩紅熱A組鏈球菌emm基因型別變化特征分析．中華流行病學雜志，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，Fuentes 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過程在痢疾疫情預測中的應用研究．中國衛(wèi)生統(tǒng)計，2016，31（3）：395-398．

［7］中華人民共和國衛(wèi)生部．WS282-2008猩紅熱診斷標準．北京：人民衛(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ī)院門診流量管理研究．中國數字醫(yī)學，2015，10（6）：77-79．

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

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

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

［14］周雨．1997-2006年沈陽市和平區(qū)猩紅熱資料分析．預防醫(yī)學論壇，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年猩紅熱疫情流行病學分析．中國公共衛(wèi)生，2012，28（6）：826-827．

［17］楊毅，趙國浩，秦愛民．面板數據的有序聚類分析及其應用-以全球氣候變化聚類分析為例．統(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．

（責任編輯：劉 壯）

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

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

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

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

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