天氣狀況因素對隼形目猛禽遷徙數(shù)量的影響
——基于負二項混合回歸模型

林芊蔚

( 愛丁堡大學數(shù)學學院, 英國愛丁堡 EH9 3FD )

猛禽是隼形目和鸮形目鳥類的統(tǒng)稱, 多為食肉型或食腐型.其中,隼形目猛禽包含鷹、鵟、隼、雕、鷲、鷂、鶚等鳥類, 主要在白天活動. 鸮形目猛禽被稱為貓頭鷹, 主要在夜晚活動. 猛禽在自然界處于食物鏈的頂端, 捕食其他鳥類和鼠、兔、蛇等, 在維持生態(tài)系統(tǒng)平衡方面發(fā)揮著重要作用.

在每年春秋遷徙季, 部分猛禽會在繁殖地與越冬地之間來回遷徙以適應因季節(jié)性變化而改變的環(huán)境條件. 為了應對遷徙途中能量大量消耗問題, 它們采取了各種減少能量消耗的策略. 大多數(shù)種類猛禽在遷徙過程中會利用熱上升氣流( 旋轉(zhuǎn)的空氣柱) 獲得升力以減少遷徙途中的能量消耗. 在此過程中溫度、風向等影響熱上升氣流形成的天氣狀況因素則可能影響鳥類的遷徙行為[1]. 此外, 研究表明, 猛禽遷徙者有利用天氣狀況信息, 提前或推遲遷徙出發(fā)時間的能力[2].

基于上述背景, 天氣狀況條件很可能影響猛禽的遷徙行為. 一些相關(guān)研究也初步探尋了兩者之間的關(guān)系. 李重和[3]等通過對中國東部沿海地區(qū)猛禽秋季遷徙現(xiàn)象進行觀察研究后認為猛禽遷徙與其途經(jīng)地的天氣條件有密切關(guān)系, 猛禽傾向選擇在晴朗、少云、可見度高、三級以上的東北風或西北風的天氣特征中進行遷徙. 李顯達[4]等通過計算黑龍江嫩江高峰林區(qū)猛禽日環(huán)志量與各氣象因子的相關(guān)系數(shù)后發(fā)現(xiàn)地面平均溫度、日最高氣溫、相對濕度與鳥類遷徙的相關(guān)性很強. 然而, 目前關(guān)于天氣狀況對猛禽遷徙行為影響的研究多為定性研究或簡單的相關(guān)性研究. 基于此, 文章嘗試使用負二項混合回歸模型(Negative Binomial Mixed Regression Model), 量化各天氣狀況對猛禽遷徙的影響, 探究天氣狀況變量與觀測到的隼形目猛禽遷徙總數(shù)以及各類隼形目猛禽遷徙總數(shù)之間可能存在的關(guān)系.

負二項混合回歸模型是廣義線性混合回歸模型(Generalized Linear Mixed Regression Model) 的一種, 適用于對具有過度離散性、且不具有獨立性的計數(shù)數(shù)據(jù)進行建模[5-6]. 此處過度離散指的是計數(shù)變量的實際方差高于其理論方差. 由于文章研究的變量隼形目猛禽遷徙總數(shù)以及各類隼形目猛禽遷徙總數(shù)均為離散的計數(shù)變量, 且經(jīng)計算其實際方差均大于理論方差,即存在過度離散性; 此外文章的計數(shù)數(shù)據(jù)是由相同或不同的觀察員觀察得出, 因而由同一觀察員觀察得到的計數(shù)數(shù)據(jù)之間可能具有相關(guān)性, 即可能存在隨機效應( 觀察者效應). 因而文章選擇使用負二項混合回歸模型來模擬天氣狀況變量對隼形目猛禽遷徙總數(shù)以及各類隼形目猛禽遷徙總數(shù)的影響.

2 研究方法

2.1 數(shù)據(jù)來源及處理方法

文章所有數(shù)據(jù)均來源于位于美國賓夕法尼亞州阿拉根尼·富壤特懸崖(Allegheny Front Escarpment) 的鳥類觀察站Allegheny Front Hawk Watch. 該地區(qū)共有16 種隼形目猛禽, 基本情況見表1. 其中蒼鷹是留鳥, 只在食物稀缺的時期進行不定向遷徙, 其余15 種都存在有規(guī)律的遷徙現(xiàn)象, 為文章研究對象. 在這15 種隼形目猛禽中, 根據(jù)生物學分類以及猛禽體型大小, 將其分為鷲、隼、雕、鷹、鶚5 類. 其中,鷲包含美洲鷲科的黑頭美洲鷲、紅頭美洲鷲;隼包含隼科的美洲隼、灰背隼、游隼; 雕包含鷹科大型鳥類金雕、白頭海雕; 鷹包含鷹科中小型鳥類白尾鷂、紋腹鷹、雞鷹、赤肩鵟、巨翅鵟、巨翅鵟、毛腿鵟; 鶚包含鶚科的魚鷹.該地區(qū)猛禽春季遷徙季為每年的2 月中旬至5月上旬,秋季遷徙季為8月中旬至12月中旬.在春秋遷徙季, 觀察站的志愿者( 包含記錄員和觀察員) 會觀察并記錄遷徙途中經(jīng)過該觀察站的猛禽的種類與數(shù)量、天氣狀況等相關(guān)信息.在每次觀測中存在1 名記錄員和1-4 名觀察員,觀察員將觀察的結(jié)果報告給記錄員, 由記錄員進行登記. 因樣本所限, 文章所使用的數(shù)據(jù)包含了2018 年春季遷徙季(2018 年2 月13 日至2018 年5 月5 日)、2018 年秋季遷徙季(2018年8 月15 日至2018 年12 月14 號)、2019 年春季遷徙季(2019 年2 月17 日至2019 年4 月20 日) 共1319 個觀察記錄( 樣本).

表1 隼形目猛禽基本情況

因文章旨在研究天氣狀況對隼形目猛禽遷徙總數(shù)以及各類隼形目猛禽遷徙總數(shù)的影響,故而選取隼形目猛禽遷徙總數(shù)以及各類隼形目猛禽遷徙總數(shù)作為因變量( 由于所有計數(shù)變量的計數(shù)對象都是隼形目鳥類, 文章設定后續(xù)部分所有因變量的名稱都省略“隼形目”3 個字),描述天氣狀況的各個指標為自變量, 因記錄的猛禽遷徙數(shù)量可能受到時間因素和觀察者因素的影響, 因而將這兩方面的變量作為控制變量考慮在內(nèi), 因而文章初始選擇的變量包含上述4 個方面, 見表2.

文章對于原始數(shù)據(jù)的處理主要包含異常值、缺失值的處理以及變量轉(zhuǎn)換處理.

首先, 原始數(shù)據(jù)中變量濕度、云層覆蓋度、氣壓的取值中存在異常值. 因上述3 個變量均為數(shù)值型變量, 對于這些變量中的異常值, 使用這些變量的均值來進行代替.

其次, 在缺失值處理方面, 變量風向和降水狀況分別存在10 例與23 例缺失值, 因這兩個變量均為分類變量, 分別使用這兩個變量的眾數(shù)來插補其缺失值. 此外, 原始數(shù)據(jù)中存在觀察員編號缺失的情況, 此時如果記錄員編號存在, 則認為該記錄員同時也是觀察員, 即可用記錄員編號插補觀察員編號. 但若同一個樣本中記錄員編號和觀察員編號都缺失, 則刪除該樣本. 經(jīng)篩選后, 數(shù)據(jù)集存在1 268 個有效樣本.

最后, 部分原始變量需要進行轉(zhuǎn)換, 以便于后續(xù)分析. 其中, 原始變量風向共包含17 個取值： 正北、北東北、東北、東東北、東、東東南、東南、南東南、南、南西南、西南、西西南、西、西西北、西北、北西北、變化. 文章將其縮減為5 個取值： 即偏北方向、偏南方向、偏西方向、偏東方向及變化風向. 此外,變量觀察時段是基于原始數(shù)據(jù)中觀察記錄的起始時間得出, 將其劃分為6：00-8：59、9：00- 11：59、12：00 -14：59、15：00 -17：59 四個時段. 另外,變量主要觀察員編號是基于原始數(shù)據(jù)集中各觀察員編號得出. 由前述可知每個樣本中觀察員有1-4 名, 若觀察員數(shù)量只有1 名, 則該觀察員是主要觀察員, 若觀察員的數(shù)量大于1 名,則選擇在所有觀察記錄中出現(xiàn)頻率最高的觀察員作為主要觀察員.

表2 變量概要表

2.2 負二項混合回歸模型方法介紹

文章所使用的負二項混合回歸模型的基本原理如下： 在非混合的負二項回歸模型中, 計數(shù)變量( 因變量) 遵循負二項式分布.

其中, yi代表第i 個計數(shù)變量或因變量, ui代表第i 個計數(shù)變量的期望,θ 代表離散參數(shù),Γ(·) 代表伽馬函數(shù), n 代表樣本數(shù). 在負二項回歸模型中, 計數(shù)變量的期望和方差如下：

從上述期望和方差的表達式中可以看出在負二項回歸模型中, 計數(shù)變量的方差總是大于期望, 離散參數(shù)θ 決定了過度離散的程度, 當θ 趨近于正無窮時, 計數(shù)變量的期望近似于其方差, 模型趨近于泊松分布, 此時不存在過度離散情況.

負二項回歸模型的前提假設條件要求因變量之間是相互獨立的, 若該假設條件不滿足,則模型需要進行改進, 其中一種方法是引入隨機效應, 用于解釋因變量之間的相關(guān)性. 引入隨機效應后( 即模型轉(zhuǎn)變?yōu)樨摱椈旌匣貧w模型), 模型的連結(jié)函數(shù)g(·) 如下：

其中, β 為固定效應向量, γ 為隨機效應向量, Xi和Zi分別為固定效應和隨機效應的設計矩陣中的對應行. 此外, 通常假定遵循多元正態(tài)分布, 其中D 是相應的方差矩陣.

此外, 計數(shù)變量可能取決于決定事件發(fā)生機會數(shù)量的規(guī)模變量(size variable), 例如觀察不同城市的入室盜竊數(shù)量, 則該數(shù)量將取決于這些城市的住戶數(shù)量, 城市的住戶數(shù)量即為規(guī)模變量[7-8]. 此時可用計數(shù)變量除以規(guī)模變量的比率作為模型的因變量進行擬合, 則該模型的連結(jié)函數(shù)可等價轉(zhuǎn)化：

2.3 模型比較方法

在模型擬合過程中, 不僅需要確定模型的類型, 還需要使用模型間比較的方法對各模型進行比較和選擇. 其中, 一種常見的方法是使用貝葉斯信息準則 (BIC).BIC 統(tǒng)計量的計算公式為：

其中, α 為模型中參數(shù)的數(shù)量, n 是樣本大小, In(n) 是模型復雜度的懲罰項, L 是模型的極大似然函數(shù).BIC 準則可用于比較不同模型的表現(xiàn), 具有較低BIC 統(tǒng)計量的模型要優(yōu)于較高BIC 統(tǒng)計量的模型.

3 模型擬合過程與結(jié)果分析

3.1 因變量描述性統(tǒng)計

在模型擬合前, 先初步觀察各個因變量的分布特點, 見圖1. 由圖1 可知, 各因變量的分布均為離散分布, 且存在大量零計數(shù), 在之后的分析中需予以考慮.

3.2 模型擬合具體過程

在觀察因變量的分布特點后, 文章開始擬合模型, 思路如下： 從空模型開始逐步添加自變量, 采用極大似然法估計負二項混合回歸模型中離散參數(shù)和連結(jié)函數(shù)中各固定效應及隨機效應的參數(shù), 觀察各參數(shù)的顯著性水平, 并結(jié)合BIC 準則進行模型間的比較, 判斷各個變量( 包含具有固定效應和隨機效應的變量) 是否應被納入模型. 其中, 變量觀察總時長是規(guī)模變量, 因為觀察的時間越長, 觀察到的猛禽遷徙總數(shù)就會越多. 最終擬合的各個模型結(jié)果如表3.

在上述過程中, 當引入主要觀察員編號的隨機效應時, 總模型、鷲模型、隼模型、鷹模型的BIC 值顯著減小, 因而這些最終模型包含隨機效應. 而雕模型與鶚模型的BIC 值增加,因而這兩個最終模型不包含隨機效應.

3.3 研究結(jié)果與分析

基于模型擬合結(jié)果, 進一步將最終模型的固定效應參數(shù)估計值及相應的置信區(qū)間可視化, 如圖2 所示.

由模型的估計結(jié)果可知, 在各天氣狀況變量中, 溫度是影響最廣的變量, 除雕遷徙總數(shù)以外, 溫度對其余5 個因變量均有顯著的正向影響, 但影響程度均較低. 具體而言, 當溫度升高1℃, 其余變量保持不變時, 因變量猛禽遷徙總數(shù)、鷲遷徙總數(shù)、隼遷徙總數(shù)、鷹遷徙總數(shù)、鶚遷徙總數(shù)對數(shù)的期望值將分別增加0.023、0.043、0.058、0.028、0.058 只. 除了溫度以外, 影響最廣的變量是可見度, 對除隼遷徙總數(shù)和鶚遷徙總數(shù)以外的4 個因變量均有顯著的正向影響, 但影響程度較低. 具體來說, 當可見度提高1 km, 其余變量保持不變時, 猛禽遷徙總數(shù)、鷲遷徙總數(shù)、雕遷徙總數(shù)、鷹遷徙總數(shù)對數(shù)的期望值將分別增加0.037、0.036、0.026、0.038 只. 此外, 風向也對猛禽遷徙總數(shù)、雕遷徙總數(shù)、鷹遷徙總數(shù)產(chǎn)生影響. 以偏北風為基準水平, 偏東風和偏南風在3 個模型中均使得因變量的對數(shù)期望值顯著增加, 其中偏東風的增長幅度大于偏南風. 但偏西風僅對雕遷徙總數(shù)有正向影響, 變換風向?qū)Ω饕蜃兞坑绊懢伙@著. 最后, 風速在隼模型和雕模型中對其因變量產(chǎn)生顯著的正向影響. 當風速提高1 級, 其余變量保持不變時, 隼遷徙總數(shù)、雕遷徙總數(shù)對數(shù)的期望值將分別增加0.435、0.196 只.

續(xù)表

圖2 模型各變量固定效應參數(shù)估計

同時, 在各控制變量中, 月份和觀察時段也各自在一些模型中對其因變量產(chǎn)生影響. 例如以2 月份為基準, 其余變量保持不變, 猛禽遷徙總數(shù)、鷲遷徙總數(shù)、雕遷徙總數(shù)、鷹遷徙總數(shù)的期望或?qū)?shù)值的期望在9 月份均有顯著增加. 又例如以時段6：00-8：59 為基準, 其余變量保持不變, 猛禽遷徙總數(shù)、雕遷徙總數(shù)、鷹遷徙總數(shù)在時段9：00-11：59、12：00-14：59 的期望值也均有顯著增加.

進一步, 將4 個含有隨機效應的模型中的隨機效應可視化, 如圖3 所示.

由圖3 可知, 各主要觀察員對于4 個不同模型因變量的觀察存在偏高或偏低的傾向. 例如當53 號觀察員為主要觀察員時, 則其對猛禽遷徙總數(shù)、鷲遷徙總數(shù)、鷹遷徙總數(shù)的觀察有偏高的傾向, 對隼遷徙總數(shù)的觀察有偏低的傾向. 具體來說, 若53 號觀察員為主要觀察員,則猛禽遷徙總數(shù)、鷲遷徙總數(shù)、鷹遷徙總數(shù)對數(shù)的期望值將比總體平均水平的提高0.400、1.212、0.514 只, 隼遷徙總數(shù)對數(shù)的期望值將比總體平均水平的減少-0.251 只. 此外, 當7號、12 號、14 號、34 號、36 號、40 號、49 號、52 號觀察員為主要觀察員時, 4 個模型的因變量都存在被偏低觀察的傾向, 而當28 號觀察員為主要觀察員時, 4 個模型的因變量都存在被偏高觀察的傾向.

此外, 可進一步計算4 個模型的組內(nèi)相關(guān)系數(shù)(Intraclass Correlation Coefficient, ICC), 結(jié)果見表3. 該統(tǒng)計量用于測量各同一類型( 此處指同一主要觀察員) 的因變量之間的相似程度, 即可分別計算同一主要觀察員觀察到的4 個因變量自身的相關(guān)性. 組內(nèi)相關(guān)系數(shù)的取值為0-1,約接近1 表示相似程度越高. 由表3 可知, 總模型、鷲模型、鷹模型的組內(nèi)相關(guān)系數(shù)分別為0.065,0.279,0.17,0.091. 其中鷲模型的組內(nèi)相關(guān)系數(shù)最高(0.279), 但也相對較小, 表明雖然隨機效應( 主要觀察員效應) 確實存在, 但同一主要觀察員觀察到的鷲遷徙總數(shù)之間的相關(guān)性水平較低. 總體而言, 雕模型和鶚模型不存在主要觀察員效應, 其余4 個模型均存在主要觀察員效應, 但該效應對各計數(shù)變量的影響均較低.

知識獲取體現(xiàn)受審核企業(yè)從認證機構(gòu)汲取知識的能力，順承著知識轉(zhuǎn)移，是認證機構(gòu)和受審核企業(yè)ISO14001實施的紐帶。認證機構(gòu)和受審核企業(yè)間知識存在異質(zhì)性和互補性，企業(yè)對審核員稀缺性知識資源的需求為ISO14001的實施打下了基礎(chǔ)。知識獲取的成果主要取決于受審核企業(yè)的學習能力及落實程度。學習能力與現(xiàn)有的知識積累有關(guān)，成正比關(guān)系。

圖3 模型隨機效應參數(shù)估計

3.4 模型檢驗

模型擬合后, 對于最終確定的模型, 其存在的合理性需要進行檢驗. 在接下來的分析中,將分別進行各個模型的多重共線性檢驗、零計數(shù)模擬檢驗以及模型擬合效果檢驗.

3.4.1 多重共線性檢驗

負二項混合回歸模型的假設包含自變量之間的獨立性, 其理想情況是各個自變量之間不存在多重共線性( 模型中的自變量間存在高度線性相關(guān)) 現(xiàn)象. 通常, 人們使用方差膨脹系數(shù)(VIF) 來衡量多重共線性水平, 若變量的VIF 值越高, 則表明它與其他變量間具有越高的相關(guān)性. 一般而言, 方差膨脹系數(shù)＜ 5 可認為自變量間多重共線性程度低.

計算各模型( 除鶚模型, 因該模型只有一個自變量, 無多重共線性) 各個自變量的方差膨脹系數(shù), 如圖4 所示, 各自變量的方差膨脹系數(shù)均＜ 5, 說明各自變量之間不存在明顯的多重共線性, 模型通過檢驗.

3.4.2 零計數(shù)模擬檢驗

前文圖1 顯示各個因變量均存在大量零計數(shù). 在此情況下需要對最終擬合的模型是否可以很好地反應大量零計數(shù)的特點進行檢驗. 若最終模型無法很好地反應該特點, 例如最終擬合的模型中零計數(shù)顯著地小于實際的零計數(shù),則可能需要考慮零膨脹模型, 該模型用于處理數(shù)據(jù)中存在過量零值的情況[9].

文章的做法如下： 將各個最終模型估計的參數(shù)值作為負二項分布的參數(shù), 隨機生成滿足該負二項分布的100 個與原數(shù)據(jù)具有相同樣本數(shù)的數(shù)據(jù)集, 計算各數(shù)據(jù)集中零計數(shù)的個數(shù),并繪制零計數(shù)個數(shù)的直方圖, 然后分別計算所有新生成數(shù)據(jù)集的平均零計數(shù)個數(shù)和實際數(shù)據(jù)的零計數(shù)個數(shù), 如圖5 所示.

由圖5 可知, 各個最終模型的平均模擬零計數(shù)個數(shù)與數(shù)據(jù)集的實際零計數(shù)個數(shù)均非常相近, 表明各最終模型已經(jīng)很好地反應了數(shù)據(jù)集中存在大量零計數(shù)的特點, 因此無需引入新模型例如零膨脹模型來處理多余的零計數(shù).

圖4 模型各自變量多重共線性檢驗

圖5 模型零計數(shù)個數(shù)分布

3.4.3 模型擬合效果檢驗

圖6 顯示各個最終模型的大部分觀測值均被較好地擬合, 即實際值與模型的擬合值非常接近, 各個模型整體的擬合效果較好.

圖6 模型擬合值與實際值比較

4 研究結(jié)論與局限性

4.1 研究結(jié)論

根據(jù)前文的分析, 部分天氣狀況變量( 溫度、可見度、風向、風速) 顯著影響隼形目猛禽遷徙總數(shù)及各類隼形目猛禽遷徙總數(shù), 但影響程度均較低, 且影響的顯著性、影響的程度在不同模型存在差別. 同時, 實證結(jié)果表明另一部分天氣狀況變量(濕度、氣壓、云層覆蓋度、降水情況) 對各因變量并沒有顯著影響. 此外,主要觀察員效應在多數(shù)模型中存在, 表明不同的觀察員會對觀察結(jié)果( 即觀察到的隼形目猛禽遷徙總數(shù)以及各類隼形目猛禽遷徙總數(shù)) 產(chǎn)生影響, 但影響程度均較低.

4.2 局限性

文章成功應用了負二項混合回歸模型擬合具有過度分散的計數(shù)數(shù)據(jù)隼形目猛禽遷徙總數(shù)及各類隼形目猛禽遷徙總數(shù), 并且檢驗主要觀察者效應的存在性. 但研究過程可能存在以下不足.

4.2.1 文章所使用的15 種隼形目猛禽數(shù)據(jù)中,12 種的遷徙類型均是部分遷徙( 即同一種猛禽中, 一部分是留鳥, 一部分是候鳥或旅鳥). 因該觀察站所處的地理位置為美國東部猛禽遷徙的咽喉要道, 因而文章假設途經(jīng)該觀察站的猛禽均是正在遷徙的. 但事實上, 由于部分留鳥的存在, 可能途經(jīng)該觀察站的猛禽并不全是正在遷徙的, 這部分數(shù)量的猛禽無法排除,因而可能會對估計結(jié)果產(chǎn)生影響, 造成一定程度的偏差.

4.2.2 在實際情況中, 天氣狀況變量之間可能是相互聯(lián)系的, 因而雖然在實證檢驗中另一部分天氣狀況變量( 濕度、氣壓、云層覆蓋度、降水情況) 對各個因變量的影響并不顯著, 但事實上它們可能會對猛禽遷徙總數(shù)產(chǎn)生一定影響, 這可能需要進一步研究.