手足口病是由肠道病毒引起的儿童常见传染病之一，并主要威胁5岁及以下儿童^[1]。近4年以来，浙江省手足口病发病率和死亡率均位居39种法定传染病前列，给儿童的健康造成了严重危害。研究表明，手足口病发病存在较为明显的季节波动和病原转换规律^{[2, 3]}，表现为肠道病毒优势血清型随着季节变动而发生转换；同时，由肠道病毒71型(EV71)所造成的手足口病具有每2～3年出现一次流行的现象^{[4, 5, 6]}，而EV71感染所导致重症和死亡的概率高于柯萨奇病毒A组16型(Cox A16)和其他肠道病毒^[3]。因此，如能以合适方法对手足口病病病原演变规律进行拟合，同时预测未来一段时间手足口病的优势病毒血清型，则将对疫情规模估计和防控措施的开展提供帮助。

手足口病病原学和发病、死亡、重症监测结果来源于《中国疾病预防控制信息系统》，资料收集统计时间段为2010年1月至2014年12月，其中各项指标均按月份进行统计。资料中，2010年1月至2014年10月数据为训练样本，2014年11月和12月数据为预测测试样本。

聚类分析又称集群分析，是依照事物的数值特征，对观察对象进行归类，其中系统聚类法是实际工作中使用最多的一种方法，其基本思想是：首先定义样品间距离及类与类之间的距离，开始时将每个样品各看成一类，将距离最近的两类合并，然后重新计算新类与其他类的距离，将距离最近的两类合并，直至所有的样品都合并成一类为止^{[7, 8]}。本研究采用系统聚类法中的类平均法对手足口病病原构成进行聚类分析。

马尔可夫链预测模型是一种动态随机数学模型，它是基于马尔可夫链，根据事件的目前状态预测其未来时刻变动状态的一种预测方法。一个n阶马尔可夫链由n个状态集合和一组转移概率所确定。该过程的任一时刻只能处于一个状态。如果在时刻t，过程处于状态S_j，则在t+1时刻它将以概率P处于状态S_k。其预测是根据状态之间的转移概率来推测系统未来的发展变化^[9]。

转移概率的公式为：P_jk⁽ⁱ⁾=M_jk⁽ⁱ⁾/M_j，其中P_jk⁽ⁱ⁾为状态S_j经过i步转移到状态S_k的转移概率，M_jk⁽ⁱ⁾为从状态S_j经过i步转移到状态S_k的次数，M_j为状态S_j出现的次数。则状态转移概率矩阵为：

因此，状态转移概率矩阵P⁽ⁱ⁾反映了系统从某种状态转移到各种状态的可能性大小，即可通过状态转移概率矩阵来预测系统未来的发展变化。

①以月度手足口病各病原比例为聚类变量(包括EV71、Cox A16和其他肠道病毒)，用聚类分析方法对月度手足口病病原构成进行聚类，以聚类后的类别作为时间点状态划分的提示依据，并结合流行病学知识确定时间点状态，并按时间顺序建立状态变化列表；②计算转移概率并建立转移概率矩阵；③利用转移概率矩阵预测状态转移；④确定未来一段时间内最有可能出现的手足口病病原构成比例，并以此估计重症死亡的规模；⑤根据2014年11、12月手足口病病原构成、重症、死亡数的实际发生情况验证预测的准确性。

手足口病发病规模以月发病数、月重症数和月死亡数3个指标反映，病原构成比例比较用χ²检验，发病数、死亡数及重症数比较用秩和检验(Kruskal-Wallis test)，统计分析用SAS 9.1软件完成。

自2010年1月至2014年10月，共58个样本月数据参与计算，聚类变量为每个月EV71、Cox A16和其他肠道病毒所占比例，根据类平均法系统聚类图中显示的各样品间距离和类别关系，再结合流行季节手足口病病原构成比例变化以及相应的重症、死亡病例出现频率的疾病特点^{[2, 3, 4, 5, 6]}，将月度手足口病病原构成分为6类，见图 1，每一类对应的手足口病病原构成比例如表 1所示，类间差异有统计学意义(χ²=4968，P<0.001)；类别1中，以Cox A16病原为优势病原，类别2则以EV71和Cox A16共同居多，类别3和类别4均以EV71病原为主，但类别4中EV71病原比例优势更为明显，类别5和类别6则以其他肠道病毒为主，但类别6中EV71和Cox A16比例较类别5上升较多。

图 1 月度手足口病病原构成比例聚类结果(OB1~58代表 58个样本月) Fig. 1 Categories of monthly HFMD pathogen proportion(OB1~58 signify all the sample month)

图选项

不同类别手足口病死亡数(χ²=23.58，P<0.001)和重症数(χ²=23.81，P<0.001)差异有统计学意义，类别5中，死亡数和重症数比例最低，而类别3和4的重症和死亡比例则明显高于其他类别；不同类别手足口病发病数(χ²=5.7，P=0.336)差异无统计学意义。见表 2。

根据各月状态建立手足口病病原构成的一步状态转移矩阵，见表 3、4。因2014年10月手足口病病原构成为状态6，由转移矩阵可得出，2014年11月手足口病病原构成处于状态6的概率约为88.9%(8/9)，而处于状态3和5的概率均为5.6%，12月手足口病病原构成处于状态6的概率为82.6%，处于状态1~5的概率分别为0、0.8%、8.1%、3.2%和5.3%；因此根据转移矩阵预测2014年1112月手足口病病原构成分别有88.9%和82.6%的可能以其他肠道病毒为主，同时EV71和Cox A16也处于相对较低的比例，出现死亡病例概率低，而出现1～3例重症病例的可能性仍较大，约为59.7%(此处概率由Poisson概率公式计算得出)。

根据实际监测结果，2014年1112月手足口病病原构成中，其他肠道比例分别为69.5%和58.3%，均以其他肠道病毒为主，EV71比例分别为13.8%和20.9%，Cox A16的比例分别为16.7%和20.8%，均无死亡病例，重症数分别为2例和4例，由此可见，预测结果与手足口病病原构成、重症数、死亡数实际情况吻合，实际数值与状态6的均值水平接近，与其他状态均值水平差异明显。

研究指出，应用马尔可夫模型预测的结果与数据状态划分密切相关，对样本数据序列的取值范围进行状态划分会因取值上下限值不同而导致不同的转移概率矩阵，预测结果也将不一致^[9]。本研究通过聚类分析方法进行状态划分，力求最大限度地降低人工经验判断状态而导致的偏倚，从而提高预测的准确性。

由聚类分析得到的手足口病病原构成类别/状态中，各类不同病原的构成比例差异明显，其中类别5和6中，其他肠道病毒处于优势地位，而类别1~4中，EV71病原的比例逐渐升高。研究显示，EV71导致重症和死亡的概率高于Cox A16和其他肠道病毒^[3]，在本研究中，与各类别对应的手足口病死亡数和重症数均值存在的差异也体现了这一规律，以其他肠道病毒比例最高的第5类中，死亡数和重症数最少，而随着EV71比例的升高，重症数和死亡数也相应明显上升。尽管EV71的传播力强于Cox A16，其基本复制指数(R0=5.48)高于Cox A16(R0=2.50)^[10]，但因手足口病发病受气候等其他因素影响较大^{[11, 12]}，可能导致了各类别间发病数差异无统计学意义。

根据马尔可夫预测结果显示，2014年11、12月手足口病病原构成较大可能以其他肠道为主，死亡出现概率低，但仍将有较少的重症病例发生，该预测结果与2014年11、12月手足口病的实际数值相吻合，预测准确性较好。根据转移矩阵和手足口病的流行特征可推测^[2]，未来一段时间，即2014年11月至2015年2月间，浙江省手足口病病原构成将保持在以其他肠道为主的态势，但根据5年各月的状态变化表显示，23月常存在病原构成状态变换的情况，比较明显的是2011年和2014年间，均由以其他肠道为主的构成状态变换为以EV71为主的构成状态，尽管这种变换发生的概率较低(5.6%)，但若在冬春季节转换时发生病原构成向以EV71为主的情况变换，则来年春夏季将很大概率保持以EV71为主的状态(如矩阵中所示，若进入状态3，则未来一段时间保持在状态3和4的概率将达64.3%)，极有可能导致疫情的高发和重症、死亡病例的增多。

马尔可夫预测不适合长期的预测分析，预测矩阵也并非固定不变，实际工作中应根据监测结果的 变换对矩阵进行调整，再进行下一阶段的预测分析，同时也结合疫情的总体特征和病原变化情况进行综合判断，从而得到准确性更高的预测结果，为疫情防控提供科学依据。