基于二值型數據的榮成天鵝湖菲律賓蛤仔空間分布與環境因子關系研究?

董建宇,胡成業,王學鋒,楊曉龍,張秀梅??

(1. 廣東海洋大學水產學院,廣東 湛江 524088; 2. 浙江海洋大學水產學院,浙江 舟山 316022)

物種的時空分布格局與環境因子之間的關系及其相互作用,一直都是生態學研究的熱點[1-3]。如何正確理解和定量揭示環境因子對物種空間分布的影響,尤其是海洋生物的時空分布與環境之間的互作關系,一直都是生態學家共同關注的問題[4-6]。然而,傳統的海洋生物調查數據(包括漁業調查數據),尤其針對目標物種開展的獨立調查數據,通常包含很多零值(即物種豐度為0)[7-8],難以滿足正態分布的假設,而很多數學模型方法恰恰又是建立在數據滿足正態分布假設的前提之上。

廣義線性模型(Generalized linear model,GLM)和廣義加性模型(Generalized additive model,GAM)等方法出現,為定量分析物種時空分布與環境因子之間的關系提供了技術手段[9]。GLM和GAM模型擴展了普通線性模型的一般假設,不再拘泥于目標物種數據必須滿足正態分布的假設,并且能夠更加靈活地處理物種與環境因子之間的關系,因此其越來越多地應用于海洋生物尤其是魚類與棲息環境的關系研究當中[10-12]。事實上大多數相關研究主要集中于魚類的資源豐度或單位捕撈努力量漁獲量(Catch per unit effort,CPUE)以及漁場分布與時空、環境因子的關系之上[13-15],鮮有關于底棲貝類的相關研究[16]。此外,雖然有學者探索利用其他類型數據結構如二項分布、負二項分布或Tweedie分布等揭示物種的時空分布與棲息環境之間的關系[17-19],但是卻鮮有利用二值型的簡化資源豐度數據(即服從二項分布)來探究棲息環境對底棲貝類空間分布影響的研究報道[20]。

菲律賓蛤仔(Ruditapesphilippinarum)營埋棲生活,廣泛分布于中國南北海區,是中國四大海水養殖貝類之一[21]。作為重要的海洋經濟貝類物種,雖然對于其生物學等方面的研究較多,但是卻鮮有關于其棲息環境條件與空間分布關系方面的研究報道。本研究將簡化的菲律賓蛤仔資源空間分布二值型數據(0表示物種豐度等于0,1表示物種豐度大于0)作為響應變量引入到GLM和GAM模型之中,首先對比分析了兩種模型對菲律賓蛤仔空間分布的擬合效果與預測性能;其次,擇優選擇GAM模型探究了山東榮成天鵝湖菲律賓蛤仔空間分布與環境因子之間的關系,以及各因子對其空間分布概率的影響。以期為深入了解天鵝湖菲律賓蛤仔分布動態和漁業管理提供參考,同時也為加強二值型漁業調查數據的應用提供科學依據。

1 材料與方法

1.1 研究區域與調查方法

天鵝湖又名月湖,位于山東半島的最東端,是一個典型的海岸潟湖,具體見圖1。天鵝湖南、北、西三面被陸地包圍,東面則由一細長條形的沙壩將其與黃海分離;東南角約80 m的潮汐汊道,是天鵝湖與外海進行水體交換的唯一途徑[22]。天鵝湖中生長著鰻草(Zosteramarina)和日本鰻草(Zosterajaponica)兩種海草,構成了天鵝湖所特有的海草床生態系統。按照經緯度將天鵝湖劃分成5″×5″的網格,每個網格對應一個采樣站位,并綜合考慮水深、海草密度和人類活動等因素將其分為A—E 5個區域。根據每個區域的面積比,采用成層等比例采樣的方法,按照A∶B∶C∶D∶E=4∶4∶5∶8∶3的比例,每次調查在A—E 5個區域中按比例隨機選取站位,共計24個[23]。分別于2017年冬季(2月)、春季(5月)、夏季(8月)和秋季(11月)對天鵝湖菲律賓蛤仔的空間分布及其生境狀況開展調查。調查的環境因子包括水深(Depth)、水溫(T)、鹽度(Sal)、溶解氧(DO)、pH、葉綠素a(Chla)、總有機質含量(TOM)。水溫、鹽度和溶解氧及pH采用YSI(YSI—650,美國)測定。葉綠素a和總有機質的測定參照《海洋調查規范》(GB/T 12763.6—2007)執行[24]。由于天鵝湖的水深受潮汐漲落影響較大,因此各調查站位的水深采用趙鵬等[25]根據衛星遙感數據反演的結果。使用開口面積0.05 m2的抓斗式采泥器采集菲律賓蛤仔,每個站位重復采樣4次。

圖1 山東半島榮成天鵝湖菲律賓蛤仔研究區域Fig.1 Study area of Manila clam Ruditapes philippinarum in Swan Lake,Rongcheng,Shandong Peninsula

1.2 模型方法

1.2.1 廣義線性模型 廣義線性模型(GLM)是對普通線性模型的推廣,在GLM模型中,響應變量可以是滿足指數分布族中的任一分布,而不再僅局限于正態分布;GLM模型通過連接函數的方式建立起響應變量的條件均值與解釋變量之間的函數關系[26]。GLM模型的一般表達式為

g(μY)=β0+β1X1+…+βiXi+ε。

(1)

式中:g(μY)為連接函數,μY為響應變量Y的條件均值;Xi為解釋變量;β0為截距,βi為回歸系數;ε為殘差。

1.2.2 廣義加性模型 廣義加性模型(GAM)是GLM模型的半參數擴展(Semi-parametric extension),該模型的唯一假設為函數是可加的且是光滑的[26]。與GLM模型類似,GAM模型使用連接函數建立起響應變量的條件均值與解釋變量的平滑函數之間的關系。GAM模型的一般表達式為

g(μY)=α+f1(X1)+…+fi(Xi)+ε。

(2)

式中:g(μY)為連接函數,μY為響應變量Y的條件均值;fi(Xi)為Xi的樣條平滑函數,Xi為解釋變量;α為截距;ε為殘差。

1.3 模型評估

采用受試者工作特征曲線(Receiver operating characteristic curve,ROC)及其下方的面積(Area under ROC curve,AUC)對模型的擬合精度進行評估[27],ROC曲線愈接近左上角其下方所圍成面積越大,即AUC值就越大,模型的精確度就越高。AUC取值范圍為0～1,AUC值越接近于1表明模型精確度越高;AUC≤0.5表明模型精確度與隨機判定結果無異;0.50.9,表明模型精確度優秀;AUC=1,表明模型精度完美[23]。

本研究中以天鵝湖各采樣站位菲律賓蛤仔豐度的簡化二值型數據作為響應變量Y(1代表資源豐度大于0,0代表資源豐度等于0);以logit函數為連接函數進行GLM和GAM模型擬合,誤差分布均為二項分布。將數據隨機分成訓練集和測試集:80%的數據用于模型訓練評估,20%的數據用于模型測試驗證。采用方差膨脹因子(Variance inflation factor,VIF)對預測變量進行共線性檢驗,剔除VIF大于3的變量[28]。所有通過VIF檢驗的變量代入模型中,進行變量篩選,以卡方(χ2)檢驗各變量的顯著性(P<0.05),以赤池信息準則(Akaike information criterion,AIC)衡量模型的擬合優度[29]。綜合比較所得GLM和GAM模型的評估驗證表現及其預測能力,選擇二者中表現較好的模型,探究菲律賓蛤仔空間分布與各項因子之間的關系。

數據分析和模型的構建過程均在R語言(R 4.2.1)[30]中完成;其中GLM和GAM模型均在“mgcv”包中實現[26];模型擬合精度和預測能力的評估通過“pROC”包[31]完成。

2 結果

2.1 菲律賓蛤仔的資源分布

天鵝湖四個季度菲律賓蛤仔的資源豐度(log2(x+1)轉化)如圖2所示。因為按照分層隨機取樣的方式進行調查,所以每個季度的調查站位不同。總體而言,天鵝湖菲律賓蛤仔的資源豐度呈現出由西北向東南方向明顯增多的趨勢,且資源豐度主要集中在位于天鵝湖東南部區域,其他區域資源豐度分布較少。

(豐度經過log2(x+1)轉化。Abundance transformed by log2(x+1))圖2 天鵝湖四個季度菲律賓蛤的豐度分布Fig.2 Abundance distribution of Ruditapes philippinarum in Swan Lake in four seasons

2.2 最優模型

以80%的數據對模型進行訓練,7個因子(水深、水溫、鹽度、溶解氧、pH、葉綠素a和總有機質)經過篩選最終分別各有3個因子入選GLM和GAM模型;其中入選GLM模型的因子為水深、溶解氧和鹽度,入選GAM模型的因子為水深、葉綠素a和總有機質含量,具體如表1所示。對于GLM而言,水深對菲律賓蛤仔的空間分布影響最大,其次為溶解氧和鹽度;對GAM而言,總有機質含量對菲律賓蛤仔的空間分布影響最大,其次為水深和葉綠素a。

表1 GLM和GAM最佳擬合模型參數Table 1 Parameters of the optimal model for GLM and GAM

GLM模型的最優擬合形式:Y～Depth+DO+Sal;GLM模型的偏差解釋率為44.8%,AIC值為64.40。

GAM模型的最優擬合形式:Y～s(Depth)+s(Chla)+s(TOM);其中s()為樣條平滑函數;GAM模型的偏差解釋率為49.9%,AIC值為64.45。

2.3 模型評估驗證與預測性能

GLM和GAM訓練模型評估與驗證結果如圖3和表2所示。實線表示訓練模型的ROC評估曲線,其對應的虛線表示模型驗證的ROC曲線;曲線下方的面積表示對應于曲線的AUC值。GAM模型的評估與驗證ROC曲線均位于GLM模型對應曲線的上方,且GAM訓練模型評估與驗證的AUC值均大于0.9,模型精確度評價等級為優秀;而GLM訓練模型評估與驗證的AUC值分別為0.89和0.77,模型精度評價等級為一般。上述結果表明,所構建的GLM和GAM模型均滿足適用條件,而GAM模型的精確度要優于GLM模型。

表2 GLM和GAM模型評估與驗證結果Table 2 GLM and GAM model evaluation and validation results

圖3 GLM和GAM模型的ROC曲線評估結果Fig.3 ROC curve evaluation results of GLM and GAM models

在GLM和GAM模型中,二值型響應變量Y的條件均值μY的取值范圍為[0,1],以最適μY值(閾值)為臨界點,模型對于響應變量Y的擬合精度最高;由表2可知,GAM模型在閾值μY處的敏感性與特異性均大于GLM模型。

以最適閾值μY對隨機選取的20%調查站位(n=19)進行預測,結果顯示:GLM模型與GAM模型對于其中13個豐度為0的站位預測結果相同,預測準確率均為100%;然而對于另外6個豐度大于0的站位,兩個模型預測結果存在較大差異,GAM模型成功預測了4個,而GLM模型僅預測到了1個,這表明GAM模型的預測性能優于GLM模型。

2.4 各因子對菲律賓蛤仔空間分布的影響

由于GAM模型在偏差解釋率、評估驗證和預測性能等方面均優于GLM模型,因此本研究選取GAM模型探究天鵝湖菲律賓蛤仔的空間分布與各因子之間的關系。水深、總有機質含量和水體葉綠素a含量均對天鵝湖菲律賓蛤仔的空間分布有顯著影響。其中總有機質對菲律賓蛤仔空間分布的單獨影響最大,水深和葉綠素a次之。各因子對天鵝湖菲律賓蛤仔空間分布的影響效應如圖4所示。天鵝湖菲律賓蛤仔空間分布與水深、總有機質含量呈負相關關系,隨著水深的增加和總有機質的升高,菲律賓蛤仔的空間分布逐漸下降;天鵝湖菲律賓蛤仔空間分布與水體葉綠素a濃度呈正相關關系,當葉綠素a濃度大于2 mg·m-3時,其對菲律賓蛤仔空間分布有明顯的正效應。

(f(x)為平滑函數,陰影表示95%置信區間。f(x) is smooth function,shaded represent 95% confidence intervals.)圖4 環境因子對菲律賓蛤仔空間分布的影響Fig.4 Effects of environmental factors on the spatial distribution of Ruditapes philippinarum

2.5 各因子對菲律賓蛤仔分布概率的影響

在GAM模型中,對于二值型響應變量Y,模型輸出的各因子平滑函數是建立在連結函數logit轉換尺度之上的。因此,為了確定各因子對菲律賓蛤仔空間分布概率的影響,需要對連接函數進行逆變換。經過逆變換后的各因子對菲律賓蛤仔資源空間分布概率的影響結果如圖5所示,隨著水深增加,菲律賓蛤仔的分布概率逐漸降低,當水深超過1.5 m后,菲律賓蛤仔的分布概率迅速減小,這與實地調查中發現的菲律賓蛤仔主要分布于天鵝湖中部偏東的淺水區域的結果相一致。隨著總有機質含量的增加,菲律賓蛤仔分布概率的變化情況可大致分為三個階段:當總有機質含量小于3.5%時,其分布概率保持在0.8以上的平穩階段;當總有機質含量大于3.5%而小于4.5%時,分布概率處于急速下降階段;當總有機質含量超過4.5%時,分布概率漸進于0的極低階段。隨著水體中葉綠素a含量的增加,菲律賓蛤仔的分布概率不斷增大,增大速率先快后慢,逐漸趨于平緩。

圖5 菲律賓蛤仔出現概率與環境因子的關系Fig.5 Relationship between the occurrence probability of Ruditapes philippinarum and environmental factors

3 討論

3.1 GLM和GAM模型

GLM和GAM模型擴展了一般線性模型的假設,具備處理更多種分布類型生態數據的能力,且能夠與線性建模和方差分析完美結合,因此,自問世以來就被廣泛地應用于生態學研究當中[9,11]。但是在解決實際生態學問題時,通常會面臨數據難以服從正態分布,且難以確定響應變量與解釋變量之間的關系,而不恰當的假設往往會導致后續數據分析的偏差,得出錯誤的結論。漁業調查數據尤其是針對目標種類開展的獨立調查往往在很多站位并未觀測到該目標種的出現,從而導致大量零值的存在[7-8]。為了滿足模型假設和減少大量零值的影響,多數研究的通常做法是將數據進行轉換使其服從正態分布,然后進行建模[32-34];然而,卻很少有關于響應變量數據服從其他類型分布如Gamma分布、負二項分布或二項分布以及Tweedie分布等研究報道[17-19]。二值型漁業調查數據是對CPUE或資源豐度數據的一種簡化,直接以1和0表示其有和無,而不考慮其量的多少,雖然會損失原始CPUE或資源豐度數據所包含的部分信息,但是數據的獲取卻相對更容易。本研究結果也表明簡化的二值型數據同樣可以很好的解釋環境因子對于物種空間分布與出現概率的影響。

GLM模型是建立在響應變量的均值與解釋變量的線性組合之間的假設關系之上的[26],但是在分析具體問題時,并非每一種解釋變量都是線性的,例如本研究中總有機質含量與響應變量Y。因此這可能是導致本研究中GLM模型解釋率相對較低和預測能力相對較差的原因。GAM模型能夠通過平滑函數較好地實現對非線性因子的擬合,訓練模型的偏差解釋率高于GLM模型,對于響應變量Y的預測準確性也從GLM模型的73.68%提高到了89.47%,這也凸顯了GAM模型在處理非線性關系中的優勢[26,35]。因此,越來越多的生態學家在解釋物種與環境之間關系時,更傾向于使用由數據所驅動的GAM模型[36-37]。

3.2 菲律賓蛤仔空間分布與環境因子之間的關系

物種的空間分布并不是完全由某單一因子所決定的,而是多個因子共同作用的結果。最優GLM模型選擇的因子偏重于水體理化參數(如溶解氧和鹽度)對菲律賓蛤仔的影響,而最優GAM模型選擇的因子偏重于沉積物理化參數即總有機質和食物需求(以葉綠素a表征)對菲律賓蛤仔空間分布的影響。入選最優GLM模型和GAM模型的因子不同,也從側面反映了菲律賓蛤仔的空間分布受多種環境因子的影響。而不同因子解釋率的高低反映了其對菲律賓蛤仔空間分布的影響程度,即解釋率越高,影響程度越大,解釋率越低,影響程度越小。水深在GLM和GAM模型均表現出較高的偏差解釋率,表明其對天鵝湖菲律賓蛤仔的空間分布的影響較大。天鵝湖中水深同海草分布的關系密切,水深較深的區域通常海草生長較為密集,形成連續的海草床。Tsai等[38]研究表明,日本鰻草的生長會降低生活于其中的菲律賓蛤仔的生存環境條件,導致菲律賓蛤仔的生長速度降低[39],此外菲律賓蛤仔的埋棲深度在日本鰻草的生長區域內也普遍較淺,更貼近于底表[38]。天鵝湖中的海草主要分布于西部區域且以鰻草占絕對優勢。鰻草的根系相較于日本鰻草更發達,因此對菲律賓蛤仔潛沙埋棲的阻礙作用也更強。此外,在天鵝湖西部海草分布區域,沉積物中的重金屬含量也相對更高[40],較高的重金屬含量會對菲律賓蛤仔產生毒害作用,因此也會限制菲律賓蛤仔的分布。

沉積物總有機質含量對天鵝湖菲律賓蛤仔空間分布的解釋率最高,表明其對菲律賓蛤仔的相對影響最大,且這種影響是負影響。盡管菲律賓蛤仔也攝食部分表層沉積物再懸浮顆粒有機物[41],但是過高的總有機質含量會導致沉積物有機質富集,對菲律賓蛤仔的生長和存活產生不利影響。沉積物總有機質含量與底質類型相關,沙含量較高的底質中總有機質含量通常較低[42]。隨著總有機質含量增加,菲律賓蛤仔的出現概率降低,這一結果恰好與其喜棲于含沙量相對較高的底質中的習性相吻合[21,43]。菲律賓蛤仔屬于濾食性雙殼貝類,主要以浮游植物為食。葉綠素a含量在一定程度上反映了水體中浮游植物的豐富度,隨著葉綠素a含量的增加菲律賓蛤仔的出現概率增大,這與其對于食物的需求相一致。物種自身的運動能力也會影響物種的空間分布。菲律賓蛤仔在其早期生活史階段營浮游生活,經過半個月左右才轉入到底棲階段并營埋棲生活[21],鑒于其移動能力較弱,因此在本研究中未考慮其遷移行為對空間分布的影響,認為其被采集到的區域即為其終年生活之所。周年實地調查均未在天鵝湖西部區域采集到菲律賓蛤仔樣本,這暗示該區域并不適合菲律賓蛤仔生活。這與先前對菲律賓蛤仔在天鵝湖的潛在適宜生境的研究結果相一致即其適生生境范圍主要集中于天鵝湖的東南部區域[23]。此外,物種的空間分布除了受非生物因子的影響外,在很大程度上也受制于不同物種之間的相互作用[19,36],今后的研究需要進一步探究種間相互作用(如捕食關系)對物種空間分布的影響。

4 結語

本研究根據山東榮成天鵝湖4個季度菲律賓蛤仔的資源豐度二值型簡化數據和環境因子數據,應用GLM和GAM模型首次探究了天鵝湖菲律賓蛤仔的空間分布與環境因子的關系,結果顯示:菲律賓蛤仔的空間分布受多種環境因子的顯著影響,其空間分布概率與葉綠素a呈正相關關系與水深、總有機質含量呈負相關關系;GAM模型在處理二值型數據方面優于GLM模型,能更好地揭示天鵝湖菲律賓蛤仔空間分布與環境因子的關系。研究結果為有效利用簡化的漁業數據提供了案例支持。