王燕芬,王旭峰,赵继军
(青岛大学复杂性科学研究所,山东 青岛,266071)
图1 2009年~2016年手足口病发病数时间序列图Fig.1 Time series of the incidence of HFMD from 2009 to 2016
本文使用的手足口病监测数据为2009年至2016年月报告发病数,来源于公共卫生科学数据中心(www.phsciencedata.cn),全国各省年人口总数及年出生率来源于中国统计局(www.stats.gov.cn)。气候数据包括月平均气温、降雨量和相对湿度,来源于中国气象信息中心(data.cma.cn)。
寒暑假时间的选取参照各省教育厅公布的历年中小学开学及放假时间通知,1月20日至2月20日为寒假时间,7月、8月为暑假。春运共计40天,一般为春节前15天和节后25天,由于每年时间不同,根据最早开始时间和最晚结束时间,选择1月中旬至3月中旬作为春运期间。
本文首先建立时间序列SIR (Time Series Susceptible Infected Recovered,TSIR)模型估算各省手足口病月传染率,模型中的参数由马尔科夫蒙特卡洛方法(Markov chain Monte Carlo,MCMC)进行估计,然后利用线性回归模型分析月平均气温、降雨量、相对湿度、开学放假和春运对各省手足口病传染率的影响。
(1)
其中,B(t)为出生人数,Yr(t)为报告发病人数,ρ为报告率,β(t)为传染率,本文假设每年同一个月的传染率相同。由于经疾病监测系统上报的发病人数仅占实际发病数的一定比例,在计算下一时刻易感者人数时候,需要除去实际感染者人数,所以需要计算相应报告率来纠正发病人数。本文使用平滑样条法对各地报告率进行估算,模型中自由度选取经验值2.5[23,27]。
实际生活中,疾病的发病人数通常服从泊松分布或负二项分布[19,25],本文假设TSIR模型中感染者人数服从泊松分布:
Y(t)~Possion(λ)
(2)
其中,λ为感染力,且λ的期望值为:
E(λ)=β(t)Yα(t)X(t)
(3)
其中,α(0<α≤1)为人群混合程度,当α<1,说明人群是不完全混合的,即部分人之间的接触高于他们与其他人群的接触;当α=1,表示人群是均匀混合。因为TSIR模型中步长的选择(两周或是一个月)不影响最终结果[13],而我们收集的数据为月发病数,所以采用时间步长为1个月。利用TSIR模型,本文估计了各省、自治区月手足口病传染率β(t)、人群混合程度α、易感者人数X(t)及报告率ρ共15个参数。
TSIR模型中的参数由马尔科夫蒙特卡洛方法(Markov chain Monte Carlo,MCMC)进行估计[28],马尔科夫链用Gibbs采样法产生,各地区初始易感者人数X(0)、人群混合程度α和传染率β(t)的先验参数都服从均匀分布,初始易感者人数X(0)取值范围为总人口的1%~40%,人群混合程度α取值范围为[0.5, 0.999]。为了保证各省相应的基本再生数在1至100范围内,其传染率初值的取值范围不同。在初始迭代1 000次后,马尔科夫链迭代1 000万次,再对最后5 000个有效参数进行采样,为了避免采样结果自相关,采样间隔取50。由于手足口病疫苗于2016年才投入使用,且对发病数影响并不明显,所以TSIR模型中可以不考虑免疫措施的影响。
最后,在获得各省传染率的基础上,估算各省手足口病基本再生数,利用线性回归模型分析月平均气温、降雨量、相对湿度、开学放假和春运期间对各省传染率的影响。气候因素为月数据变量,开学放假(“学期”、“寒假”或“暑假”)和春运(“春运期间”或“平常时期”)均为分类变量。由于估计的传染率为月数值,所以本文选用2月作为寒假时间,7月、8月作为暑假时间。每年春运时间并不固定,考虑到人群接触后疾病的传播存在一定的时间滞后,本文选取有两个星期滞后的2月和3月作为线性回归模型中的春运期间。我们首先只选取气候因素进行线性回归,为了避免各因素之间相关性导致的模型中多重共线性问题,对变量进行筛选,将方差膨胀因子(方差膨胀因子越大,共线性越严重)大于5的因素依次去掉后,再找出对手足口病具有显著影响的气候因素。其次,应用包括具有显著影响的气候因素、开学放假和春运期间在内的线性回归模型,分析各因素对手足口病传染率的影响。
利用马尔科夫蒙特卡洛方法估算TSIR模型中的参数,对最后5 000个有效参数进行采样(图2),传染率β(t)、人群混合度α和易感者人数X(t)收敛在一定区间内。
图2 以湖南省为例,MCMC-TSIR模型最后5 000次迭代采样结果Fig.2 Taking Hunan Province as an example, the results of the last 5 000 iterations of the MCMC-TSIR model
中国各省手足口病传染率的峰值有明显的季节性(见图3),且在2月都有较大的增幅,这可能与春运期间大规模人口流动引起的人群接触率快速增加有关,因为没有发现气候因素在2月有所变化。在暑假期间,多数省的传染率处于较低水平。为了更好进行描述,根据传染率峰值发生时间,可以将全国各省归类到为4个区域(见表1):
图3 全国各省手足口病传染率峰值分布及各区域传染率偏移值Fig.3 Peak distribution of HFMD transmission rate and the deviance of β(t) in each region
地区包含省份东南地区海南、广东、广西、浙江、江苏、上海、安徽、云南、四川、重庆、贵州、山西、陕西、北京、天津、山东、江西、湖北、湖南、福建、河北、河南北部地区内蒙古、吉林、辽宁和沈阳西北地区新疆、青海、甘肃和宁夏西藏西藏
属于东南地区的各省(除山东省和北京市外)的手足口病传染率峰值集中于2月~3月;属于西北地区的省份传染率峰值在4月;东南地区和西北地区各省的手足口病传染率在秋季还有一个小高峰。北方地区传染率峰值在5月;西藏传染率峰值则在8月,这可能是因为西藏地区手足口病的平均感染年龄在5岁左右,即学龄儿童更易感染疾病,因此学校开学、放假导致的学生间接触率的变化能够影响其传染率的季节性。
根据模型估计的结果(见图4):全国多数省份手足口病基本再生数在15至35之间,只有广东省的手足口病基本再生数高于40,贵州省基本再生数则小于10;手足口病报告率在各地区也有较大差异,从整体上看,处于东南地区的省份手足口病报告率高于北部、西北地区,其中山西省、青海省、新疆和西藏自治区的报告率不足5%,广西省和海南省的报告率却高达35%,其余省份的报告率多处于10%至25%之间;各省手足口病易感者比例则集中在10%至35%之间。
图4 MCMC-TSIR模型估算各省参数结果Fig.4 Results of estimated parameters in different province of MCMC-TSIR model
利用MCMC-TSIR模型所得参数估计发病数,对实际报告发病数和预测发病数取对数后进行分析。模型预测的发病数和实际发病数之间的关系为图5。
手足口病传染率季节性在不同地区受不同因素影响(见表2):东南地区各省份传染率季节性仅与春运有关;北部地区手足足病传染率季节性则受平均气温(解释贡献率20.10%)、暑假(解释贡献率49.42%)和春运(解释贡献率8.11%)的影响;西北地区传染率季节性主要受到相对湿度(解释贡献率31.27%)和春运(解释贡献率9.80%)的影响;西藏省手足口病传染率季节性既不受气候因素影响,也不受春运的影响。
图5 模型估计发病数及报告发病数关系Fig.5 The relationship between the estimated number of cases and the number of reported cases
影响因素估计值平均温度相对湿度暑假春运期间R2东南地区回归系数---1.204e-060.642 2P值-< 2.2e-16北部地区回归系数 7.565e-08--1.944e-069.748e-070.420 3P值4.93e-06-0.000 1630.075 522西北地区回归系数--1.141e-07 -1.710e-060.414 9P值-0.010 85-0.068 26西藏回归系数-----P值----
各地区手足口病均不受降雨量影响,因此表中去除了降雨量这一项。地区内各省存在差异性,对各省的手足口病传染率季节性影响因素进行分析发现东南地区各省传染率季节性都与春运有显著相关,其中广东省传染率季节性还与相对湿度有关,北京、天津传染率季节性还与相对湿度和温度相关;北部地区各省的传染率季节性都与相对湿度有关,吉林、辽宁和黑龙江省传染率季节性还受平均温度影响;属于西北地区的宁夏、甘肃手足口病传染率季节性与相对湿度、温度有关,青海省传染率仅与春运有关,新疆自治区传染率季节性则与相对湿度和降雨量有关。
手足口病传染率在中国各省都有明显的季节性,且在2月都有大幅度增加,这与春运期间,整体人群接触率变化有关。根据传染率峰值的时间分布,可以将全国各省归类为四个区域:东南地区、西北地区、北部地区及西藏。东南地区的多数省份传染率高峰集中在2月和3月,西北地区,包括新疆、青海、甘肃和宁夏,手足口病传染率高峰发生在4月,北方地区,包括内蒙古、辽宁、吉林和沈阳的传染率在5月达到最大值,而西藏地区的传染率在8月达到峰值。
根据MCMC-TSIR模型,全国多数省份手足口病基本再生数在15至35之间,这可能与人口密度、易感者比例等因素有关。根据估算的基本再生数,手足口病的免疫覆盖率应高于93.3%。从整体上看,多数省份的报告率在10%至25%之间,处于东南地区的省份手足口病报告率高于北部地区和西北地区,而西藏自治区报告率仅2.5%,这可能与当地医疗水平、公共卫生健康监测系统的完善程度有关;各省手足口病易感者比例则集中在10%至35%之间。
手足口病传染率的季节性影响因素主要分为两类:社会因素及气候因素。东南地区的大多数省份的手足口病传染率季节性仅受春运期间的显著影响。东南地区多数省份经济发达,流动人口数目高于其他地区,春运期间,人口迁移量更是急剧增加,使得各年龄组的接触率远高于平常时间,封闭、拥挤的旅途环境使得疾病的传播可能性大大增加。北方地区的传染率季节性则与平均温度和暑假相关,西北地区的传染率季节性受到相对湿度和春运的影响;这可能是因为西北、北部地区的月气候变化较为显著,气候变化会影响病毒在环境中的生存率,还会在一定程度上影响人的行为方式,间接影响人群接触率变化。西藏地区传染率不受气候、学期和春运的影响。在寒暑假期间,各省的手足口病传染率也有所增加,尤其西藏地区,其手足口病传染率峰值在8月。这可能是因为中国各省手足口病平均感染年龄不同,西藏地区手足口病的平均感染年龄在5岁左右,学龄儿童更易感染疾病[4],因此学校开学、放假导致的儿童间接触率的变化对其传染率的季节性的影响较其他省更加显著。
为了更好了解手足口病传染率的季节性及其影响因素,需要更进一步研究,在未来研究中,我们将具体分析不同年龄组手足口病传染率季节性及其影响因素,并将影响因素代入模型中,通过曲线下的面积(Area Under Curve,AUC)来进行评估。