气质量达标（API＜100）天数和月平均相对湿度以及月平均气温5个变量与人感染H7N9禽流感之间的相关度，利用针对地理人口因素的空间相关性分析和对气象因子的时间维度上的非线性分析，量化各危险因子对疾病发生的影响，建立人感染H7N9禽流感预警模型。 资料与方法

1. 数据来源：全国人感染H7N9禽流感病例的发病日期和暴露地点资料（截至2014年6月1日）源自国家卫生和计划生育委员会网站^[8]、WHO网站^[9]和Flutrackers论坛^[10]。地市级年均人口数量、行政区域面积、空气质量达标天数等数据来自中国城市统计年鉴^[11]，农贸市场数量来源于高德地图（www.amap.com），气象资料（2013年2月至2014年1月）包括月平均气温和月平均相对湿度来自中国气象科学数据共享服务网（http：//cdc.cma.gov.cn）。

2. 建立模型：本研究采用空间自回归（SAR）模型^{[12, 13, 14]}。并基于最小二乘回归（OLS）基础上的空间相关性检验，选用空间滞后模型进行分析。其表达式

y＝ρWy＋x₁ β₁＋x₂ β₂＋x₃ β₃＋ε

式中y为各地级市发病的空间平滑率（spatial rates，以邻近地区为单位的发病率计算）；x₁、x₂、x₃代表活禽市场密度、人口密度和空气质量达标天数3个危险因子； ρ为空间回归系数； Wy是空间滞后因变量；β₁、β₂、β₃为待估参数；ε为随机误差。本研究的空间分析基于Queen相邻型空间权重矩阵^[15]。 采用广义相加模型（GAM）分析气温和相对空气湿度非线性变量^[16]。本研究假定第t月人感染H7N9禽流感日报告病例数Yt 近似服从Poisson分布，期望值为μt 。其连接函数为对数函数，则表达式

Yt～Poisson（ μt）

ln（μ_t）＝f1（x_1t）＋f₂（x_2t）＋f₃（seas）

式中Yt是某地级市第t月的人感染H7N9禽流感报告病例数；f1（x_1t）、f2（x_2t）、f₃（seas）是自然样条函数（natural spline）^[17]，其中x1t、x2t分别是对应地级市第t月平均气温和平均相对空气湿度，seas是基于一年中月份的序列数据，以控制模型的季节性^[18]，同时使各地级市的发病例数据保持时间上的可比性^[19]。采用GeoDa 1.6.5、 R 3.1.1软件计算。

3. 模型整合与危险预测：SAR模型得出的各地区发病率预测值反映该地区长期存在的发病风险，而GAM得出的预测值则反映各地区整体发病的动态趋势。一个地区特定时间的预测发病风险由2个模型的预测值共同组成（图 1）。

图 1 我国各地级市人感染H7N9发病风险的双模型计算

图选项

将预测的月平均气温及其相对湿度（基于上一年同期数值）代入GAM，可得到危险预测值，并与SAR模型的风险预测值分别标准化后，共同计算最终的地区发病风险。预测模型的建立是基于前75%的病例样本量（330例），并采用后25%的发病数据（110例）作为预警模型的验证。 结 果

1. 发病概况：截至2014年6月1日我国人感染H7N9禽流感疫情除8例定位不详病例，本研究总病例样本量为440例。用于模型建立的330例（前75%的病例）其空间分布于55个地级市，主要集中在我国东部沿海地区（图 2）；发病时间呈现双峰分布，即2013年3、4月和2014年1月有明显的发病高峰（图 3）。

图 2 2013－2014年我国330例人感染H7N9禽流感病例的空间分布

图选项

图 3 2013－2014年我国330例人感染H7N9禽流感病例的发病时间分布

图选项

2. 危险因子计算： （1）SAR模型：对活禽市场密度、人口密度和空气质量达标天数3个危险因子利用OLS得到其估计结果（表 1），其中活禽市场密度（P＜0.001）与人口密度（P＝0.008）对发病的影响有统计学意义。空间自相关检验结果（表 2）显示我国各地级市之间人感染H7N9禽流感发病率存在很强的空间自相关（Moran’s I＝22.42，P＜0.001）。空间滞后模型估计结果（表 3）显示，新加入的空间回归系数相关性显著（P＜0.001）。对模型进行检验，空间滞后残差值的Moran’s I统计量在检验水平下无正负相关关系（Moran’s I＝0.04，P＝0.038 0），表明空间滞后模型已经排除了空间自相关的情况，对数据的空间相关分布拟合良好。

表 1 采用OLS模型估计我国人感染H7N9禽流感的相关危险因子

表选项

表 2 OLS模型的空间自相关诊断结果

表选项

表 3 利用空间滞后模型估计我国人感染H7N9禽流感的相关危险因子

表选项

（2）GAM：本文气温和相对空气湿度2个变量的平滑自由度选择为4，季节性平滑自由度选择为3，建立的模型能较合理解释潜在风险而不出现过度拟合。模型的变量曲线拟合结果（表 4）显示各变量的拟合显著。图 4显示各变量的曲线拟合情况。月均气温变量的拟合曲线（图 4?）显示气温5～20 ℃时曲线拟合值维持在高位；月均相对空气湿度与拟合值之间的曲线关系（图 4?）显示高危值集中在50～65；发病的季节性较为明显，即1－12月逐月下降且呈近似线性（图 4?）；为更直观体现2个变量对病例预测值的共同影响，本研究构建了对病例预测值的三维图形（图 4?）。

表 4 GAM的曲线拟合结果

表选项

图 4 GAM拟合曲线图

图选项

3. 危险地区预测模型的构建和验证：发病危险预测（综合以上2个模型）按发病风险的标准差描绘成标准差地图，实际发病情况（后25%病例）按发病例数的具体值描绘成独特值地图（图 5）。对比危险预测地图与实际发病地图发现，2月发病的26个地级市全部位于地图的红色危险区域内，预测成功率为100%；3月的模型预测成功率为86.7%（15个实际发病城市中有13个位于发病预测危险区域中）；4月模型预测成功率为66.7%（18个实际发病城市中有12个位于发病预测危险区域中）；5月模型预测成功率为63.6%（11个实际发病城市中有7个位于发病预测危险区域中）。

图 5 2014年2－5月我国人感染H7N9禽流感危险预测和实际发病情况

图选项

人感染H7N9禽流感疫情具有复杂的因果链。从染疫禽类输入、染疫禽类间病毒扩散，到人群接触染疫动物、暴露人群感染发病，各环节中均存在多个危险因素。此外，人感染H7N9禽流感还受到各地不同疾病监测水平、实验室检测水平和健康事业投入等社会因素的影响。本研究选取活禽市场密度、人口密度2个危险因子主要衡量人群暴露的机会，气温和相对湿度主要评价病毒扩散的风险，并与空气质量共同作用于暴露人群的感染和发病。研究中未考量信息偏倚等社会因素的影响。考虑到我国不同地区发生人感染H7N9禽流感疫情存在明显的空间及季节性聚集，本研究将上述发病相关危险因子分别拟合进入2个模型，评估空间和时间与发病的相关性，再将两模型结合作出最终的危险预测模型。

在OLS模型中活禽市场密度、人口密度两变量与空间平滑发病率表现出显著的相关，回归系数均为正值。并与Fang等^[20]的研究结果一致，表明活禽市场密度或人口密度较高的地区，人感染H7N9禽流感发病风险较高，对相关地区应加强疫情防控。空气质量达标天数与发病情况之间的相关关系并不明显，可能是该项指标不能精确反映当地实际的空气质量水平，或是“空气质量”这一变量对发病的直接作用相对较弱。在空间滞后模型中空间回归系数有显著相关性，模型拟合度较OLS提高，说明中国人感染H7N9禽流感疫情具有很强的空间相关性。但是其他回归系数的显著性均出现绝对值的下降，表明 “活禽市场密度”、“人口密度”对发病的作用在某种程度上是通过空间邻近关系而表达。 在寒冷（5～20 ℃）和相对空气湿度低（50%～65%）的条件下，人感染H7N9禽流感发病风险较高，与以往研究结果相吻合^{[20, 21]}。季节性变量显示在一年中的前6个月是发病高危月份，可见对于我国大部分地区而言，冬春季的气温和相对空气湿度是发病的高危险因素，此期间应更加警惕人感染H7N9禽流感疫情的暴发。

本研究将两模型得出的危险预测值标准化后相加，获得地区的特定发病风险。标准化的过程消除了两模型预测值之间量纲的差异，将两模型的预测值相加意味着两者的权重分配均等。在2014年2月的危险预测地图中3个高发病风险地带与实际发病聚集地（月病例数＞2例）对应，验证了高发病风险区域。济宁、菏泽、六安、云浮4个地级市标记为低发病风险区域，但与其接壤的所有其他地级市则被标记为高发病风险区域。实际情况证实无发病的4个地级市周围出现了不同程度的发病地区，表明该模型可在短期内准确预测所有的低发病风险区域。随着时间的推移，模型预测的准确度逐月下降，但仍能准确预测发病趋势。模型预测我国人感染H7N9禽流感发病高风险地区逐步从南方向北方移动，而实际上2014年4、5月，在吉林省延边和山东省烟台地区分别发生1例和2例病例，验证了模型预测的准确性。

模型中北京地区有4个月一直被标记为高发病危险预测区域，但实际上始终未出现病例。笔者认为这是由于北京地区的孤立发病情况在SAR模型中赋予了较高的空间权重，而积极的控制措施使疫情未出现扩散，故该地区其后的发病例数始终为零。与此相似，在危险预测模型的后3个月中，我国西部地区由于处于适宜的发病温度和湿度，GAM对这些地区的预测值增加，标记为高发病危险区域，然而实际上由于人群暴露风险较低，即使气象条件符合，也不易发生疫情。

综上所述，本研究拟合出的模型能在短期内较准确地预测人感染H7N9禽流感发病地区，具有实用意义。而建立更精确的危险预测模型还有赖于纳入更多的发病影响因素，以及选择更高的时空分辨率和更准确的方法。