基于二值型数据的荣成天鹅湖菲律宾蛤仔空间分布与环境因子关系研究❋

董建宇,胡成业,王学锋,杨晓龙,张秀梅❋❋

(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模型,能更好地揭示天鹅湖菲律宾蛤仔空间分布与环境因子的关系。研究结果为有效利用简化的渔业数据提供了案例支持。