基于ZINB模型和气象因素的禽霍乱发病预测

高 翔，肖建华，王洪斌

（东北农业大学黑龙江省动物普通疾病防治重点实验室， 哈尔滨 150030）

0 引 言

在养殖业中，家禽养殖受传染性疾病影响较大。与发达国家相比，中国家禽养殖在疾病防控方面存在的一个问题是，从业者往往更依赖疫苗与抗生素进行传染病的防控，忽视养殖管理、综合防制的重要性。经验证明，忽视日常管理，仅依靠抗生素与疫苗无法完全有效地控制家禽传染性疾病的传播流行[1]。疫苗虽然可以阻止动物发病，却不能阻止在养殖场与养殖家禽种群内已存在的病毒的扩散与传播[2]。同时，由于病原微生物不断的变异，在面对新发传染病时，疫苗的研发会有一个滞后期。在这种情况下，评估养殖场内外存在的影响疫病传播的风险因素，针对性加强管理，提高动物抗病能力，降低动物接触传播病原的风险，是一种有效的防控措施。同时，抗生素残留产生的耐药问题已严重影响对相关疾病的治疗。抗生素随着动物排泄物进入环境，在土壤和水体中慢慢富集，也最终对人类健康造成危害。对中国河流湖泊水样进行抽样调查后发现，在采集到的水体样本中检出的抗生素浓度普遍较高[3]。

综上所述，国内外流行病学研究专家越来越重视利用现代信息技术加强动物养殖健康管理[4-5]。同时，传染性疾病的传播流行往往受到环境地理等因素，如气象、植被、海拔等的影响[6]。为了及早发现传染性疾病暴发流行的“苗头”，国内外的学者利用各种疾病与风险因素数据，结合数学评估预测模型，开展对动物传染性疾病风险因素评估与早期预警的研究，在实际应用中表现出了良好的性能[7-8]。

禽巴氏杆菌病(Avian pasteuretlosis)，俗称禽霍乱，是一种在家禽中传播的急性败血性呼吸系统传染病。长久以来，禽霍乱是中国《兽医公报》中收录发病数据最多的家禽传染性疾病。由于该病发病快，致死率高，严重威胁着中国家禽养殖业的健康发展，中国将此病列入家禽传染病控制二类名单[9-10]。通过对过往发病数据分析，发现禽霍乱在中国集中分布于南方高温高湿气候省份，并表现出季节性发病的特点，因此，气象因素是影响禽霍乱发病的重要风险因素。

基于以上背景，该研究提出利用数学评估模型评估气象数据对禽霍乱发病影响，利用网络地理信息（Web GIS）技术建立禽霍乱监测预警系统的思路。通过本文中的研究与系统开发，期望可以支持家禽养殖场疫病防控管理决策，降低经济损失，为数学评估模型在各类动物疫病防控工作中的应用研究提供参考。

1 系统设计

1.1 系统结构设计

系统的建设采用 B/S架构，不仅能够满足用户对传统信息管理系统的需求，还可以进行数据的输入、输出、编辑、更改、查询、计算等功能。而且，在地理信息技术的支持下，还可以进行数据的电子地图可视化操作与空间分析功能。用户在客户端填写的各类信息通过互联网络，传输到管理部门的后台数据库服务器，进行运算分析，得到的结果再通过网络返回用户。系统总体结构可以被抽象为用户层、功能层、技术层及数据层 4个层次。

图1 系统架构图Fig.1 Architecture diagram of system

1.2 系统功能设计

系统功能大体上划分为地图基本操作、疫情数据上报、疫情数据查询、疫情空间分析和风险预测 5大功能模块。具体目标包括：1）地图基本操作：包括地图的放大、缩小与平移，鹰眼视图等功能。2）疫情数据上报：疫情上报信息包括发病日期、发病地点、动物种类、动物来源、发病数量、报告人等。利用系统，使用者点击地图上紧急疫情发生的地点，弹出数据录入对话框，录入相关信息。信息录入后，在地图上生成特定标示的红色点。数据录入还包括对历史疫情信息的收集整理，包括发病数、死亡数、捕杀数等。对于历史疫情数据，还可以形成统计报表和统计图形。3）疫情数据查询：信息查询功能包括图查属性功能与属性查图功能。图查属性功能允许使用者以感兴趣的某一点为中心，通过设置查询半径，在生成的查询区范围内的所有信息都已列表的形式被显示。这方便了使用者及时掌握疫情点周围的信息；属性查图通过键入感兴趣点的信息，查找其在地图上的位置。这方便了使用者快速查找到所感兴趣的疫情点。4）疫情空间分析：包括缓冲分析与路径分析。缓冲分析为三重缓冲分析，对应满足了疫情管理中需要划分疫区、受影响区以及封锁区的要求；路径分析通过在地图上设置起点及终点，由系统自动给出两点间最短的路径信息。方便使用者在应对紧急疫情时进行物质调配和人员的调动。5）风险预测：利用得到的各风险因素与疾病发病相关的风险系数，结合数据库中的风险因素数据，计算疾病发生的风险值并生成风险地图。地图以颜色的深浅判断疾病风险的大小。

1.3 系统数据库设计

系统数据库包括属性数据库、空间数据库以及模型库。数据库支持系统实现对数据的管理、计算功能，是系统实现空间分析与预测的基础。本研究采用的是ArcGIS 10.0内置数据库来管理空间数据，以矢量地图格式，按空间、地物属性分类的依据分层存储，地图坐标为经纬度坐标。

图2 监测预警系统细节功能图Fig.2 Function of monitoring early warning system

属性数据库主要包括养殖状况数据以及历史疫情数据等。养殖情况数据包括全国各省、市、自治区每年年底家禽存栏数、出栏数、各省的土地面积等数据，为显示疫情信息等功能提供数据支持。疫情数据包括各省月家禽疫情数据，包括发病数、死亡数、捕杀数等数据，是建立评估预测模型的主要数据来源。

预测模型是整个决策支持系统的核心功能部分。模型库主要是利用已有的专家知识和专家经验，采用零膨胀负二项式回归模型(zero-inflated negative binomial，ZINB)方法建立相关预警模型，通过编程技术建立系统模型。模型可调用数据库中的数据及参数值，并将结果返回动态数据库。

2 预测模型的建立

2.1 研究数据

该研究中，2006年1月至2016年12月中国各省禽霍乱每月发病数据来源于中国兽医局网站下载的《兽医公报》。将2006年1月至2015年12月间的数据资料用于预测研究，2016年1月至2016年12月的数据资料作为模型预测效果的检验数据。养殖数据来自于由国家统计局网站下载的2006―2016年《中国统计年鉴》。气象数据来自于国家气象局信息中心。

2.2 时间与季节性分析

该研究中，首先通过时间与季节性分析排除之前月份发病对本月发病的影响，以及季节因子对疾病发病的影响，为气象因素风险作用的分析做准备。

时间序列分析方法通过判断时间上前后两点间数据的相关关系，探寻数据相关事件的发展规律。时间序列分析包括对数据随机性与平稳性的分析。

时间序列分析的一个前提是要求研究数据为非平稳性序列。检验时间序列是否为平稳性序列，目前常用的是增广 ADF检验(augmented dickey-fuller test)，简记为ADF检验[11]。ADF检验是单位根检验(unit root test)的一种，原假设为序列具有非平稳性，显著水平取值为0.01。当P＜0.01拒绝原假设。ADF检验任何一部分结果中出现P＞0.01的情况，都表明原假设成立。其计算公式如下所示：

式中μ为常数均值；β为回归系数；εt为误差项。

进行时间序列分析的另一个前提是作为研究对象的数据序列必须存在时间上的关联性，纯随机序列没有进行时间序列分析的价值[12]。理论上，纯随机序列数据的自相关系数（autocorrelation coefficients，AC）为 0。自相关系数可根据以下公式求得：

式中 n为序列观测值的数目；ρk为样本自相关系数；k(0＜k＜n)为延迟阶数；xt为t时刻序列测量值。

该研究中，首先利用ADF方法判断禽霍乱发病数据序列的平稳性。然后，利用自相关与偏自相关系数判断发病数据本身存在的自相关性。相关分析由SAS 9.2分析软件完成。

该研究通过使用季节分解模型对2006年1月至2015年12月中国家禽传染性疾病月发病数据进行分析，描述其季节分布特征与长期趋势。通过季节分解模型，可以将原始发病数据分解得到季节因子序列（seasonal adjustment factors，SAF），通过SAF序列，可以判断发病数据的季节性，当序列峰值大于 1时，发病数据具有季节性[13]。对于时间序列ARIMA（n, d, m），其季节分解模型的一般形式为：

式中e为时间序列ARIMA（n, d, m）的残差；d为时间序列由非平稳序列转变为平稳序列的d阶差分；a为时间序列ARIMA（n, d, m）的白噪音，n与m分别为时间序列自相关系数最佳的阶层和偏自相关系数最佳的阶数，

2.3 风险分析

通过对禽霍乱发病数据时空分布特点的研究，发现禽霍乱发病会受到季节因子作用的影响。同时，上月发病会影响本月疾病的发生。因此，季节因子与上月发病数作为协变量纳入模型分析。当地的家禽养殖量作为哑变量控制因养殖规模不同而对发病数据造成的影响。

为了检验气象因子相互之间是否存在共线性。研究利用 Spearman 等级相关分析方法对气象因子间的相关关系进行检验。Spearman 等级相关系数取值区间为0～1，一般认为相关系数大于 0.7时，两变量间存在共线性，不能同时纳入多变量的模型分析[14]。

对于存在共线性的 2个气象因素，研究采用采用赤池信息量信息准则（akaike information criterion，AIC）对纳入了不同变量的模型进行评估。根据相关准则，AIC的取值越小，模型的拟合预测结果越好[15-16]。

式中K是所拟合模型中参数的数量；L是对数似然值；n是观测值数目。

本文中的研究使用互相关(cross-correlation)分析方法研究气象因素对全国家禽传染性疾病的滞后作用，其计算公式为[17]：

式中 Cov(yt,Xt-k)为延迟 k阶的协方差，Var(yt)，Var(Xt-k)分别是相应变量的方差。

在传染性疾病的研究中，由于疾病空间集聚以及季节性分布的特点，发病数的方差往往明显偏离发病均数，不能满足等离散的条件。对这些类型的传染病数据进行分析，大大超过了传统计数模型的分析能力。同时，在对发病数据此类计数资料进行研究的过程中，时常还会遇到在单位观测时间段内，事件的发生次数为0的现象。传统的计数模型对此类含有大量“0”的事件数据都不能有效的进行拟合。因此，Lambert于1992年提出了ZINB模型[18]。ZINB模型是一种将“零膨胀”事件的发生过程分作了零计数过程以及泊松计数过程 2个部分的混合概率分析方法。利用零计数过程可以解释事件发展过程中出现大量“0”的原因，有效地弥补了传统计数分析模型分析含有大量“0”的事件数据不足的问题，使得分析结果更加准确[19]。

ZINB模型形式如下：

式中φ为ZI参数；α为散度参数，E(Y)=(1-φ)/λ。

2.4 模型评估

该研究釆用组内相关系数（intraclass correlation coefficient，ICC）分析方法评估模型预测发病数与实际疾病发病数之间的一致性[20]。本文采用的ICC分析方法是经Shrout和Fleiss对传统的组内相关系数分析方法理论作出改进，将其扩展到了评价重复量测数据之间相关性应用领域后的方法，组内相关系数被定义为被测量者间变异占总变异的比例[21]。改进后，ICC分析方法被逐渐应用于经济学与医学统计分析领域，组内相关系数的取值范围在0～1，一般认为ICC大于0.75时预测值与实际值拟合良好[22]。其计算公式为：

式中MS处理为真实值与预测间的均方；MS区组代表不同月间病例数的均方；MS误差为误差的均方；n被观察对象的例数；k为处理组数。

3 预测结果与预测系统实现

3.1 时间序列与季节性分析结果

由ADF检验结果可见，各检验方法均有P值大于0.01，不能拒绝无效假设，可认为序列为非平稳性时间序列。

表1 禽霍乱发病数序列单位根检验结果Table 1 Results of ADF tests of Avian pasteuretlosis sequence

从自相关系数可以看出，除 1～4阶自相关系数外，延迟5阶及以后的自相关系数均在2倍的标准误以内，没有统计学意义。同时，由表2可以看出延迟1阶的自相关系数及偏自相关系数对比其他阶系数均明显偏大，这说明了中国2006年1月至2015年12月禽霍乱病例具有延迟一阶的自相关，即当月的禽霍乱发病数受上个月病例数的影响。

表2 禽霍乱发病数序列自相关与偏自相关系数Table 2 Autocorrelation coefficients(AC) and partial autocorrelation coefficients(PAC) of Avian pasteuretlosis sequence

通过禽霍乱各月份发病数与季节因子SAF的曲线可以看出，2006年10月至2015年12月中国禽霍乱发病多出现在秋冬季节（图 3），具有明显的季节性（季节因子峰值超过1）[13]。

图3 禽霍乱各月份发病数Fig.3 Number of Avian pasteurellosis cases in every month

3.2 气象因素风险分析

该研究将气象因素的当月数据和之前 1～3月的数据作为自变量，将全国家禽传染性疾病每月发病数作为响应变量，纳入分析。模型纳入前 1个月的家禽传染性疾病月发病数用于控制疾病自相关；纳入季节因子控制禽病疫情发生的季节性；纳入各地区的禽养殖数量总数抵消养殖数量不同造成的差异。

图4 禽霍乱的季节因子分析Fig.4 Seasonal factor analysis of Avian pasteurellosis

由于Spearman等级相关分析显示，发现中国当月的降水量与月平均湿度以及月平均温度，当月的降水量与滞后1个月的月最高温度呈共线性，不能同时纳入模型。因此，禽霍乱多因素研究有2个模型：I滞后本月的平均相对湿度+本月的平均温度+本月的最高气温+本月的平均风速+上月的平均气温；II本月的总降水量+本月的最高气温+本月的平均风速+上月的平均气温。根据AIC准则，选择模型I作为最优模型。

利用回归模型分析得到气象因素对禽霍乱发病影响作用结果如表 1示。在计数过程中，本月最高气温、本月相对湿度与本月禽霍乱发生呈正相关，本月和上月的月平均气温也与其呈正相关；季节性因子与本月平均风速与禽霍乱发生呈负相关。在零计数过程中，上个月禽霍乱病例对零计数的发生概率呈负相关，即上个月的禽霍乱病例会降低本月不发生禽霍乱的概率。

表3 ZINB回归模型结果Table 3 Results of multivariate analysis by using zero-inflated negative binomial model

3.3 模型评估结果

利用组内相关系数方法对预测值与中国2016年禽霍乱月实际发病数进行对比，结果显示 ICC取值分别为0.802，说明模型预测值与实际病例数据有较好的一致性，预测值与观测值之间的R2取值为0.7853。

图5 2016年禽霍乱的月发病数预测值与观测值拟合图Fig.5 Predicted number of outbreaks versus observed number of monthly-reported outbreaks of Avia, 2016

3.4 系统实现

系统基于ArcGIS 10.0与Visual Studio.NET 2008开发平台，以 C#. NET语言为基础进行开发，数据库采用Microsoft SQL Server 2005。开发完成的系统软件可在Windows操作系统下运行，系统界面如图6所示。系统要实现的主要目标包括：①数据管理：通过疫情上报的方式对家禽历史疫情信息包括发病数、死亡数、捕杀数等进行收集；通过系统自身建立的环境数据专题库对影响家禽传染性疾病发病的气象因素进行整理；②疾病预警：利用ZINB模型分析得到的各气象风险因素与疾病发病相关的风险系数，计算疾病发生的风险值并生成风险地图；③应急管理：通过缓冲分析、路径分析等多种空间分析功能，方便使用者在应对紧急疫情时制定物质调配和人员调动等应急管理策略；④用户可随时查询任意时间段内的监测记录，并生成相对应的报告图片。系统部分功能如图7所示。

4 讨论

传统研究中，利用逻辑回归分析或相关性分析模型揭示气象因素对传染性疾病发生传播的影响作用，在进行此类研究时，需要预先假设气象因素与传染病发病之间满足线性或对数线性相关的关系。然而，在实际中，气象因素与疾病发病之间的关系往往表现非线性特点。例如，在一定的范围内，病原微生物的致病能力随环境温度的上升而增加。但当气温超过一定范围时，环境温度升高反而抑制病毒的感染能力[23]。此外，对于传染性疾病发生发展，气象因素还存在滞后性影响。例如，当气象变化超过生物体适宜的生存范围时，机体的免疫能力会受损下降。但这种变化并不是立即发生的，生物体的内环境具有一定的“稳态”，对于外界条件的变化往往表现一定的抵抗力[24]。ZINB模型是一种改进后的泊松分析方法，它不仅能够很好地解释气象因素对传染性疾病的非线性影响，结合互相关分析，对气象因子的滞后性作出评估。而且，对于疫情数据集中因为某些月份没有疾病发生而存在的大量“零数据”， ZINB模型也能有效消除其对分析结果的影响，大幅提髙预测结果的准确度。因此，ZINB模型被越来越多的应用于气象对疾病传播流行影响的分析中[25]。

ZINB模型分析显示禽霍乱发生与气温、湿度呈正相关，与风速呈负相关。气象因素与气候变化对传染病发生的影响表现在各个方面，包括影响宿主的抵抗力和易感性等[26]。一般而言，机体内环境对于外界环境的变化存在一定的抵抗力，维持着一定的“稳态”。但是当动物长期处于不适宜的外界环境中时，内环境的“稳态”会受到破坏，引起免疫力的下降。Garren等[27]于1956年的研究表明，长期处于高温、高湿的环境下会导致机体的免疫能力受到抑制。解剖后发现，长期在高温、高湿环境下饲养的鸡的胸腺、法氏囊等免疫器官发生坏死萎缩、淋巴液相较于正常的鸡明显减少，抗病能力显著下降。

除了高温、高湿的环境影响外，家禽养殖场内及其周边环境的空气质量也会造成家禽免疫力的下降。这是由于健康的动物体内或体外也会带有一些致病微生物，这些微生物随唾液或脱落的皮屑在养殖场内或养殖场周边的空气中飘荡，经过长时间的风化，与空气中的水、颗粒尘埃等进一步融合形成气溶胶。在气溶胶内的微生物往往可以在很长时间内保持感染能力[28]。同时，家禽养殖场内及周边的空气中，硫化氢与氨气的浓度往往很高。这些气体进入呼吸道，会引起呼吸道免疫系统的纤毛-黏液系统的功能受损，降低呼吸道的抗感染能力[29]。

综上所述，高气温与湿度引起禽霍乱发病数的增加可能是因为高的气温与湿度会造成动物体免疫机能的损伤与下降。由于机体对外界环境的变动具有一定的抵抗力，因此气温对禽霍乱发病数的影响也表现出一定的滞后作用。较高的风速可能在一定程度上抵抗高温、高湿的影响，同时，风速较大时可能会降低养殖场环境空气中氨气、粉尘及有害气体的浓度。因此，风速对禽霍乱的发生呈现一定的负相关作用。

5 结论

本文通过将季节性分析、时间序列分析方法和零膨胀负二项式回归模型技术相结合，对中国禽霍乱发病的分布特征进行了描述，对气象因素对中国禽霍乱发病的影响作用进行了评估。结果表明：1）禽霍乱发生与气温、湿度呈正相关，与风速呈负相关。同时，气象因素影响的滞后作用可能反映了动物机体对外界环境变动具有一定的抵抗能力；2）基于零膨胀负二项式回归模型的分析结果，利用气象数据对禽霍乱的发病作出预测，与实际发病数比较，准确率可达 80.2%，该模型的分析结果能够较好地反映相关气象因子对中国禽霍乱发病的影响作用。

家禽传染性疾病的传播过程复杂，受多种因素影响，同时，中国家禽养殖模式多种多样，养殖水平高低不一。条件所限，本文只选择对气象因素对中国禽霍乱发病的影响进行了分析。如何将经济、交通乃至候鸟迁徙等更多风险因素纳入分析，进一步提高预测的准确度，是下一步研究工作的重点。