天气状况因素对隼形目猛禽迁徙数量的影响——基于负二项混合回归模型

林芊蔚

( 爱丁堡大学数学学院, 英国爱丁堡 EH9 3FD )

猛禽是隼形目和鸮形目鸟类的统称, 多为食肉型或食腐型.其中,隼形目猛禽包含鹰、鵟、隼、雕、鹫、鹞、鹗等鸟类, 主要在白天活动. 鸮形目猛禽被称为猫头鹰, 主要在夜晚活动. 猛禽在自然界处于食物链的顶端, 捕食其他鸟类和鼠、兔、蛇等, 在维持生态系统平衡方面发挥着重要作用.

在每年春秋迁徙季, 部分猛禽会在繁殖地与越冬地之间来回迁徙以适应因季节性变化而改变的环境条件. 为了应对迁徙途中能量大量消耗问题, 它们采取了各种减少能量消耗的策略. 大多数种类猛禽在迁徙过程中会利用热上升气流( 旋转的空气柱) 获得升力以减少迁徙途中的能量消耗. 在此过程中温度、风向等影响热上升气流形成的天气状况因素则可能影响鸟类的迁徙行为[1]. 此外, 研究表明, 猛禽迁徙者有利用天气状况信息, 提前或推迟迁徙出发时间的能力[2].

基于上述背景, 天气状况条件很可能影响猛禽的迁徙行为. 一些相关研究也初步探寻了两者之间的关系. 李重和[3]等通过对中国东部沿海地区猛禽秋季迁徙现象进行观察研究后认为猛禽迁徙与其途经地的天气条件有密切关系, 猛禽倾向选择在晴朗、少云、可见度高、三级以上的东北风或西北风的天气特征中进行迁徙. 李显达[4]等通过计算黑龙江嫩江高峰林区猛禽日环志量与各气象因子的相关系数后发现地面平均温度、日最高气温、相对湿度与鸟类迁徙的相关性很强. 然而, 目前关于天气状况对猛禽迁徙行为影响的研究多为定性研究或简单的相关性研究. 基于此, 文章尝试使用负二项混合回归模型(Negative Binomial Mixed Regression Model), 量化各天气状况对猛禽迁徙的影响, 探究天气状况变量与观测到的隼形目猛禽迁徙总数以及各类隼形目猛禽迁徙总数之间可能存在的关系.

负二项混合回归模型是广义线性混合回归模型(Generalized Linear Mixed Regression Model) 的一种, 适用于对具有过度离散性、且不具有独立性的计数数据进行建模[5-6]. 此处过度离散指的是计数变量的实际方差高于其理论方差. 由于文章研究的变量隼形目猛禽迁徙总数以及各类隼形目猛禽迁徙总数均为离散的计数变量, 且经计算其实际方差均大于理论方差,即存在过度离散性; 此外文章的计数数据是由相同或不同的观察员观察得出, 因而由同一观察员观察得到的计数数据之间可能具有相关性, 即可能存在随机效应( 观察者效应). 因而文章选择使用负二项混合回归模型来模拟天气状况变量对隼形目猛禽迁徙总数以及各类隼形目猛禽迁徙总数的影响.

2 研究方法

2.1 数据来源及处理方法

文章所有数据均来源于位于美国宾夕法尼亚州阿拉根尼·富壤特悬崖(Allegheny Front Escarpment) 的鸟类观察站Allegheny Front Hawk Watch. 该地区共有16 种隼形目猛禽, 基本情况见表1. 其中苍鹰是留鸟, 只在食物稀缺的时期进行不定向迁徙, 其余15 种都存在有规律的迁徙现象, 为文章研究对象. 在这15 种隼形目猛禽中, 根据生物学分类以及猛禽体型大小, 将其分为鹫、隼、雕、鹰、鹗5 类. 其中,鹫包含美洲鹫科的黑头美洲鹫、红头美洲鹫;隼包含隼科的美洲隼、灰背隼、游隼; 雕包含鹰科大型鸟类金雕、白头海雕; 鹰包含鹰科中小型鸟类白尾鹞、纹腹鹰、鸡鹰、赤肩鵟、巨翅鵟、巨翅鵟、毛腿鵟; 鹗包含鹗科的鱼鹰.该地区猛禽春季迁徙季为每年的2 月中旬至5月上旬,秋季迁徙季为8月中旬至12月中旬.在春秋迁徙季, 观察站的志愿者( 包含记录员和观察员) 会观察并记录迁徙途中经过该观察站的猛禽的种类与数量、天气状况等相关信息.在每次观测中存在1 名记录员和1-4 名观察员,观察员将观察的结果报告给记录员, 由记录员进行登记. 因样本所限, 文章所使用的数据包含了2018 年春季迁徙季(2018 年2 月13 日至2018 年5 月5 日)、2018 年秋季迁徙季(2018年8 月15 日至2018 年12 月14 号)、2019 年春季迁徙季(2019 年2 月17 日至2019 年4 月20 日) 共1319 个观察记录( 样本).

表1 隼形目猛禽基本情况

因文章旨在研究天气状况对隼形目猛禽迁徙总数以及各类隼形目猛禽迁徙总数的影响,故而选取隼形目猛禽迁徙总数以及各类隼形目猛禽迁徙总数作为因变量( 由于所有计数变量的计数对象都是隼形目鸟类, 文章设定后续部分所有因变量的名称都省略“隼形目”3 个字),描述天气状况的各个指标为自变量, 因记录的猛禽迁徙数量可能受到时间因素和观察者因素的影响, 因而将这两方面的变量作为控制变量考虑在内, 因而文章初始选择的变量包含上述4 个方面, 见表2.

文章对于原始数据的处理主要包含异常值、缺失值的处理以及变量转换处理.

首先, 原始数据中变量湿度、云层覆盖度、气压的取值中存在异常值. 因上述3 个变量均为数值型变量, 对于这些变量中的异常值, 使用这些变量的均值来进行代替.

其次, 在缺失值处理方面, 变量风向和降水状况分别存在10 例与23 例缺失值, 因这两个变量均为分类变量, 分别使用这两个变量的众数来插补其缺失值. 此外, 原始数据中存在观察员编号缺失的情况, 此时如果记录员编号存在, 则认为该记录员同时也是观察员, 即可用记录员编号插补观察员编号. 但若同一个样本中记录员编号和观察员编号都缺失, 则删除该样本. 经筛选后, 数据集存在1 268 个有效样本.

最后, 部分原始变量需要进行转换, 以便于后续分析. 其中, 原始变量风向共包含17 个取值： 正北、北东北、东北、东东北、东、东东南、东南、南东南、南、南西南、西南、西西南、西、西西北、西北、北西北、变化. 文章将其缩减为5 个取值： 即偏北方向、偏南方向、偏西方向、偏东方向及变化风向. 此外,变量观察时段是基于原始数据中观察记录的起始时间得出, 将其划分为6：00-8：59、9：00- 11：59、12：00 -14：59、15：00 -17：59 四个时段. 另外,变量主要观察员编号是基于原始数据集中各观察员编号得出. 由前述可知每个样本中观察员有1-4 名, 若观察员数量只有1 名, 则该观察员是主要观察员, 若观察员的数量大于1 名,则选择在所有观察记录中出现频率最高的观察员作为主要观察员.

表2 变量概要表

2.2 负二项混合回归模型方法介绍

文章所使用的负二项混合回归模型的基本原理如下： 在非混合的负二项回归模型中, 计数变量( 因变量) 遵循负二项式分布.

其中, yi代表第i 个计数变量或因变量, ui代表第i 个计数变量的期望,θ 代表离散参数,Γ(·) 代表伽马函数, n 代表样本数. 在负二项回归模型中, 计数变量的期望和方差如下：

从上述期望和方差的表达式中可以看出在负二项回归模型中, 计数变量的方差总是大于期望, 离散参数θ 决定了过度离散的程度, 当θ 趋近于正无穷时, 计数变量的期望近似于其方差, 模型趋近于泊松分布, 此时不存在过度离散情况.

负二项回归模型的前提假设条件要求因变量之间是相互独立的, 若该假设条件不满足,则模型需要进行改进, 其中一种方法是引入随机效应, 用于解释因变量之间的相关性. 引入随机效应后( 即模型转变为负二项混合回归模型), 模型的连结函数g(·) 如下：

其中, β 为固定效应向量, γ 为随机效应向量, Xi和Zi分别为固定效应和随机效应的设计矩阵中的对应行. 此外, 通常假定遵循多元正态分布, 其中D 是相应的方差矩阵.

此外, 计数变量可能取决于决定事件发生机会数量的规模变量(size variable), 例如观察不同城市的入室盗窃数量, 则该数量将取决于这些城市的住户数量, 城市的住户数量即为规模变量[7-8]. 此时可用计数变量除以规模变量的比率作为模型的因变量进行拟合, 则该模型的连结函数可等价转化：

2.3 模型比较方法

在模型拟合过程中, 不仅需要确定模型的类型, 还需要使用模型间比较的方法对各模型进行比较和选择. 其中, 一种常见的方法是使用贝叶斯信息准则 (BIC).BIC 统计量的计算公式为：

其中, α 为模型中参数的数量, n 是样本大小, In(n) 是模型复杂度的惩罚项, L 是模型的极大似然函数.BIC 准则可用于比较不同模型的表现, 具有较低BIC 统计量的模型要优于较高BIC 统计量的模型.

3 模型拟合过程与结果分析

3.1 因变量描述性统计

在模型拟合前, 先初步观察各个因变量的分布特点, 见图1. 由图1 可知, 各因变量的分布均为离散分布, 且存在大量零计数, 在之后的分析中需予以考虑.

3.2 模型拟合具体过程

在观察因变量的分布特点后, 文章开始拟合模型, 思路如下： 从空模型开始逐步添加自变量, 采用极大似然法估计负二项混合回归模型中离散参数和连结函数中各固定效应及随机效应的参数, 观察各参数的显著性水平, 并结合BIC 准则进行模型间的比较, 判断各个变量( 包含具有固定效应和随机效应的变量) 是否应被纳入模型. 其中, 变量观察总时长是规模变量, 因为观察的时间越长, 观察到的猛禽迁徙总数就会越多. 最终拟合的各个模型结果如表3.

在上述过程中, 当引入主要观察员编号的随机效应时, 总模型、鹫模型、隼模型、鹰模型的BIC 值显著减小, 因而这些最终模型包含随机效应. 而雕模型与鹗模型的BIC 值增加,因而这两个最终模型不包含随机效应.

3.3 研究结果与分析

基于模型拟合结果, 进一步将最终模型的固定效应参数估计值及相应的置信区间可视化, 如图2 所示.

由模型的估计结果可知, 在各天气状况变量中, 温度是影响最广的变量, 除雕迁徙总数以外, 温度对其余5 个因变量均有显著的正向影响, 但影响程度均较低. 具体而言, 当温度升高1℃, 其余变量保持不变时, 因变量猛禽迁徙总数、鹫迁徙总数、隼迁徙总数、鹰迁徙总数、鹗迁徙总数对数的期望值将分别增加0.023、0.043、0.058、0.028、0.058 只. 除了温度以外, 影响最广的变量是可见度, 对除隼迁徙总数和鹗迁徙总数以外的4 个因变量均有显著的正向影响, 但影响程度较低. 具体来说, 当可见度提高1 km, 其余变量保持不变时, 猛禽迁徙总数、鹫迁徙总数、雕迁徙总数、鹰迁徙总数对数的期望值将分别增加0.037、0.036、0.026、0.038 只. 此外, 风向也对猛禽迁徙总数、雕迁徙总数、鹰迁徙总数产生影响. 以偏北风为基准水平, 偏东风和偏南风在3 个模型中均使得因变量的对数期望值显著增加, 其中偏东风的增长幅度大于偏南风. 但偏西风仅对雕迁徙总数有正向影响, 变换风向对各因变量影响均不显著. 最后, 风速在隼模型和雕模型中对其因变量产生显著的正向影响. 当风速提高1 级, 其余变量保持不变时, 隼迁徙总数、雕迁徙总数对数的期望值将分别增加0.435、0.196 只.

续表

图2 模型各变量固定效应参数估计

同时, 在各控制变量中, 月份和观察时段也各自在一些模型中对其因变量产生影响. 例如以2 月份为基准, 其余变量保持不变, 猛禽迁徙总数、鹫迁徙总数、雕迁徙总数、鹰迁徙总数的期望或对数值的期望在9 月份均有显著增加. 又例如以时段6：00-8：59 为基准, 其余变量保持不变, 猛禽迁徙总数、雕迁徙总数、鹰迁徙总数在时段9：00-11：59、12：00-14：59 的期望值也均有显著增加.

进一步, 将4 个含有随机效应的模型中的随机效应可视化, 如图3 所示.

由图3 可知, 各主要观察员对于4 个不同模型因变量的观察存在偏高或偏低的倾向. 例如当53 号观察员为主要观察员时, 则其对猛禽迁徙总数、鹫迁徙总数、鹰迁徙总数的观察有偏高的倾向, 对隼迁徙总数的观察有偏低的倾向. 具体来说, 若53 号观察员为主要观察员,则猛禽迁徙总数、鹫迁徙总数、鹰迁徙总数对数的期望值将比总体平均水平的提高0.400、1.212、0.514 只, 隼迁徙总数对数的期望值将比总体平均水平的减少-0.251 只. 此外, 当7号、12 号、14 号、34 号、36 号、40 号、49 号、52 号观察员为主要观察员时, 4 个模型的因变量都存在被偏低观察的倾向, 而当28 号观察员为主要观察员时, 4 个模型的因变量都存在被偏高观察的倾向.

此外, 可进一步计算4 个模型的组内相关系数(Intraclass Correlation Coefficient, ICC), 结果见表3. 该统计量用于测量各同一类型( 此处指同一主要观察员) 的因变量之间的相似程度, 即可分别计算同一主要观察员观察到的4 个因变量自身的相关性. 组内相关系数的取值为0-1,约接近1 表示相似程度越高. 由表3 可知, 总模型、鹫模型、鹰模型的组内相关系数分别为0.065,0.279,0.17,0.091. 其中鹫模型的组内相关系数最高(0.279), 但也相对较小, 表明虽然随机效应( 主要观察员效应) 确实存在, 但同一主要观察员观察到的鹫迁徙总数之间的相关性水平较低. 总体而言, 雕模型和鹗模型不存在主要观察员效应, 其余4 个模型均存在主要观察员效应, 但该效应对各计数变量的影响均较低.

知识获取体现受审核企业从认证机构汲取知识的能力，顺承着知识转移，是认证机构和受审核企业ISO14001实施的纽带。认证机构和受审核企业间知识存在异质性和互补性，企业对审核员稀缺性知识资源的需求为ISO14001的实施打下了基础。知识获取的成果主要取决于受审核企业的学习能力及落实程度。学习能力与现有的知识积累有关，成正比关系。

图3 模型随机效应参数估计

3.4 模型检验

模型拟合后, 对于最终确定的模型, 其存在的合理性需要进行检验. 在接下来的分析中,将分别进行各个模型的多重共线性检验、零计数模拟检验以及模型拟合效果检验.

3.4.1 多重共线性检验

负二项混合回归模型的假设包含自变量之间的独立性, 其理想情况是各个自变量之间不存在多重共线性( 模型中的自变量间存在高度线性相关) 现象. 通常, 人们使用方差膨胀系数(VIF) 来衡量多重共线性水平, 若变量的VIF 值越高, 则表明它与其他变量间具有越高的相关性. 一般而言, 方差膨胀系数＜ 5 可认为自变量间多重共线性程度低.

计算各模型( 除鹗模型, 因该模型只有一个自变量, 无多重共线性) 各个自变量的方差膨胀系数, 如图4 所示, 各自变量的方差膨胀系数均＜ 5, 说明各自变量之间不存在明显的多重共线性, 模型通过检验.

3.4.2 零计数模拟检验

前文图1 显示各个因变量均存在大量零计数. 在此情况下需要对最终拟合的模型是否可以很好地反应大量零计数的特点进行检验. 若最终模型无法很好地反应该特点, 例如最终拟合的模型中零计数显著地小于实际的零计数,则可能需要考虑零膨胀模型, 该模型用于处理数据中存在过量零值的情况[9].

文章的做法如下： 将各个最终模型估计的参数值作为负二项分布的参数, 随机生成满足该负二项分布的100 个与原数据具有相同样本数的数据集, 计算各数据集中零计数的个数,并绘制零计数个数的直方图, 然后分别计算所有新生成数据集的平均零计数个数和实际数据的零计数个数, 如图5 所示.

由图5 可知, 各个最终模型的平均模拟零计数个数与数据集的实际零计数个数均非常相近, 表明各最终模型已经很好地反应了数据集中存在大量零计数的特点, 因此无需引入新模型例如零膨胀模型来处理多余的零计数.

图4 模型各自变量多重共线性检验

图5 模型零计数个数分布

3.4.3 模型拟合效果检验

图6 显示各个最终模型的大部分观测值均被较好地拟合, 即实际值与模型的拟合值非常接近, 各个模型整体的拟合效果较好.

图6 模型拟合值与实际值比较

4 研究结论与局限性

4.1 研究结论

根据前文的分析, 部分天气状况变量( 温度、可见度、风向、风速) 显著影响隼形目猛禽迁徙总数及各类隼形目猛禽迁徙总数, 但影响程度均较低, 且影响的显著性、影响的程度在不同模型存在差别. 同时, 实证结果表明另一部分天气状况变量(湿度、气压、云层覆盖度、降水情况) 对各因变量并没有显著影响. 此外,主要观察员效应在多数模型中存在, 表明不同的观察员会对观察结果( 即观察到的隼形目猛禽迁徙总数以及各类隼形目猛禽迁徙总数) 产生影响, 但影响程度均较低.

4.2 局限性

文章成功应用了负二项混合回归模型拟合具有过度分散的计数数据隼形目猛禽迁徙总数及各类隼形目猛禽迁徙总数, 并且检验主要观察者效应的存在性. 但研究过程可能存在以下不足.

4.2.1 文章所使用的15 种隼形目猛禽数据中,12 种的迁徙类型均是部分迁徙( 即同一种猛禽中, 一部分是留鸟, 一部分是候鸟或旅鸟). 因该观察站所处的地理位置为美国东部猛禽迁徙的咽喉要道, 因而文章假设途经该观察站的猛禽均是正在迁徙的. 但事实上, 由于部分留鸟的存在, 可能途经该观察站的猛禽并不全是正在迁徙的, 这部分数量的猛禽无法排除,因而可能会对估计结果产生影响, 造成一定程度的偏差.

4.2.2 在实际情况中, 天气状况变量之间可能是相互联系的, 因而虽然在实证检验中另一部分天气状况变量( 湿度、气压、云层覆盖度、降水情况) 对各个因变量的影响并不显著, 但事实上它们可能会对猛禽迁徙总数产生一定影响, 这可能需要进一步研究.