★ FEATURED ARTICLE 做生态数据分析这几年被问到最多的问题之一就是功能多样性指数到底怎么算网上资料不少但要么只讲概念不讲代码要么直接甩一段代码但数据格式对不上跑起来全是报错。尤其FRic、FEve、FDiv这几个常用指数很多人第一次用FD包就被数据格式卡住了矩阵转置、性状标准化、物种名匹配这些前置工作看着不起眼实际占了整个分析流程的一大半时间。这篇东西我尽量把坑都填上从数据整理到指数计算再到结果解读给出一套可以直接照着跑的完整流程。先说清楚这篇内容适合谁。如果你是做植物群落、动物群落、微生物群落或者任何涉及物种功能性状和多样性测度的研究想用R语言计算功能丰富度、功能均匀度、功能离散度这类指数这篇文章可以直接作为操作手册用。零基础能用有基础也能从中抠一些细节出来。我用的核心工具包是FD这是生态学里做功能多样性最主流的包另一个是vegan主要用于数据整理和前期检查。1. 功能多样性指数的生态学意义与核心概念1.1 为什么要计算功能多样性传统多样性测度比如Shannon-Wiener指数、Simpson指数看的是物种数和物种相对多度说白了就是有多少种、每种多少数量。但这类指数有个天然盲区它假设所有物种在生态系统里扮演的角色是等价的可实际上两片林地即使物种数一样多也可能一片全是阔叶树种、另一片全是针叶树种它们在养分循环、林冠结构、动物栖息地供给上的生态功能完全不同。这就是功能多样性要回答的问题——不只看你有多少物种更看这些物种的性状组合有多丰富、在功能空间里铺得有多开。举个例子一片只有单一优势种的农田和一片物种数相同但性状差异巨大的次生林物种多样性可能接近但功能多样性差异会非常大。而功能多样性恰恰与生态系统过程密切相关比如生产力、养分保持能力、抵抗入侵能力。现在的群落生态学研究尤其是环境梯度分析和人为干扰响应研究里功能多样性几乎已经是标配指标了如果审稿人看到你的文稿通篇只有物种多样性没有功能多样性大概率会提意见。1.2 FD包能算哪些功能多样性指数Laliberté和Shipley开发的FD包提供了围绕功能空间和树状图两大类框架的一系列指数。最常被顶刊引用的有以下几个指数缩写全称中文俗称生态学含义FRicFunctional Richness功能丰富度群落占据功能空间的大小凸包体积值越大说明功能策略越多样FEveFunctional Evenness功能均匀度多度在功能空间分布的均匀程度反映资源利用的均衡性FDivFunctional Divergence功能离散度多度在功能空间边缘的分散程度高值意味着极端性状的物种多度高FDisFunctional Dispersion功能分散度物种到重心加权距离的平均值群落在功能空间的整体离散程度RaoQRaos Quadratic EntropyRao二次熵两两物种间的性状距离加权值同时考虑多度与性状差异这些指数的计算基础是功能距离也就是物种间基于多个性状计算出的距离矩阵。FD包默认用欧氏距离也可以自己传入距离矩阵。距离矩阵的质量直接决定功能多样性指数的可信度所以性状选择和处理特别重要后面我会详细讲。2. 数据格式整理三条铁律与实战内功2.1 群落数据矩阵的标准格式与常见丑格式FD包做多样性计算的函数是dbFD()它要求两个核心输入x是群落数据矩阵a是性状数据矩阵。这两个对象的格式有严格的潜规则很多人的报错根源就在这儿。群落矩阵x必须满足这几个条件行是样方/群落列是物种内容是物种的某种数量指标多度、盖度、生物量等推荐用多度或盖度行名必须是样方ID列名必须是物种名且不能有重复不需要合计列或百分比列保持原始测度单位即可内部会自动标准化我见过最多的丑格式有这三种第一种物种名做了数字编码比如把样方里出现的物种编成sp1、sp2、sp3。这本身没问题问题在于编码表和性状数据的物种名对不上后面merge的时候直接全变NA。第二种群落数据是长格式即三列样方ID、物种名、多度。这种格式是野外调查的原始输出但dbFD()不认长格式必须先用tidyr::pivot_wider()或reshape2::dcast()转成宽格式样方ID做行名、物种名做列名、多度为单元格值。第三种把多个调查年份或者多个处理混在一个文件里行名就变成了2021_样方A_处理1这种组合ID。倒不是说不能算但后续如果要做分组比较就得提前拆列建议还是分开存储别在一张表里硬塞。2.2 性状数据矩阵的格式要求与标准化处理性状矩阵a的格式同样有严格要求行是物种且行名必须与群落矩阵里的物种名完全一致列是性状变量每一列必须是数值型连续变量这里重点强调每一列必须是数值型连续变量这句话。FD包不认因子型数据也不认字符型数据。如果你有一个性状是叶片质地革质/纸质/膜质这种分类变量直接放进去会报错或者在计算距离时被当数值处理结果完全错误。处理方法有两种一种是把分类性状转成虚拟变量这个建议慎重。比如叶片质地有3个水平就生成3列0/1变量但虚拟变量会显著增加性状维度小数据集里很容易把功能空间撑变形导致FRic失真。另一种更稳妥的做法是先问自己这个分类性状对研究问题真的必要吗如果它的生态功能已经被其他连续性状覆盖了比如比叶面积、叶片厚度等就直接舍弃。功能多样性计算的核心价值在连续性状的组合效应抓大放小比什么都往模型里塞更靠谱。性状数据的标准化是另一个关键点。不同性状的测量单位差异很大比如种子质量是克树高是米叶片氮含量是百分比如果不做标准化量纲大的变量会在欧氏距离里占绝对主导地位计算出来的功能距离基本被种子质量一个变量绑架了。dbFD()函数里有一个stand.x参数默认是TRUE会自动对性状矩阵做标准化。我很建议你明白它的原理而不是完全交给默认值——它做的是把每列减去均值再除以标准差即z-score标准化这样所有性状在功能空间里权重等同。2.3 名称匹配与数据清洗的实操细节名称匹配是功能多样性分析中出错率最高的一步。我复盘过很多次报错发现绝大多数问题出在以下两处第一肉眼看不见的空格。Excel里输入物种名时中文输入法或者复制粘贴经常带来首尾空格R语言里Quercus 和Quercus是两个完全不同的字符串。所以拿到数据的第一步就是用trimws()或者dplyr::mutate(across(everything(), str_trim))把空格清干净。第二同物异名和命名不统一。同一个物种群落表里写的是拉丁全名Quercus mongolica性状表里只写Q. mongolica或者群落表用的是中文名性状表用的是拉丁名。这种问题没有万能解法只能靠人工核对建议在脚本里专门加一个反查步骤算完交集的物种数和原始物种总数比一比差多少一目了然。下面的代码片段展示了一个典型的数据清洗流程我一般把它命名为数据体检每拿到一套新数据都会先跑一遍# 假设comm为群落数据矩阵trait为性状数据矩阵 library(dplyr) # 第一步检查行名列名是否干净 rownames(trait) - trimws(rownames(trait)) colnames(comm) - trimws(colnames(comm)) # 第二步找出两个数据集共有的物种 common_species - intersect(colnames(comm), rownames(trait)) common_species如果common_species的长度比实际物种数少很多就是物种名不匹配的信号必须回去核对原始表。千万别跳过这一步直接算否则看似跑出了结果实际计算只用到了共有物种那些被默默丢弃的物种里可能藏着功能上非常独特的类群。3. 动手计算功能多样性指数完整案例流程3.1 数据模拟与读取为了演示得清楚我自己模拟一个小型数据集。假设有6个样方调查到10个物种测量了4个性状比叶面积SLA、叶片氮含量Nmass、植株高度Height、种子质量SeedMass。先构造一个物种多度矩阵技术上怎么造数据都可以关键是格式要标准化set.seed(42) # 6个样方10个物种 comm - matrix( sample(0:20, 60, replace TRUE), nrow 6, dimnames list( paste0(Site, 1:6), paste0(sp, 1:10) ) ) comm注意我用sample()造出来的矩阵可能有些物种在某些样方里是0这非常贴近真实野外数据——不是每个样方都会出现所有物种。零值在群落矩阵里没问题dbFD()完全允许稀疏矩阵。再构造性状数据矩阵trait - data.frame( SLA runif(10, 10, 30), Nmass runif(10, 0.5, 3), Height runif(10, 1, 20), SeedMass runif(10, 0.1, 5) ) rownames(trait) - paste0(sp, 1:10) trait两个矩阵通过sp1到sp10完成物种对接。在实际项目中你可能还要处理多个样地的重复调查我建议把每个样方-年份对作为一个独立的行保留或者根据你的研究问题提前聚合这一步没有标准答案取决于分析设计。3.2 用dbFD()计算功能多样性指数这是核心的一步。加载FD包运行dbFD()library(FD) # 计算功能多样性指数 fdiv_out - dbFD( x comm, a trait, w rep(1, ncol(comm)) # 物种多度权重默认就是1这里显式写出来 ) # 输出结果 fdiv_out运行之后你会得到一个列表对象里面包含FRic、FEve、FDiv、FDis、RaoQ、qual.FDis等元素。qual.FDis是物种性状空间与多度分布之间的一致性度量我后面再说它怎么用。一次就能跑通是理想状态。实际报错最常见的就两类一类是提示There are zero or negative values in the a matrix——这说明性状矩阵里有0或者负值。FRic计算要取凸包体积而凸包体积的计算涉及对性状距离取对数变换FD包里有一个在内部把矩阵开根号处理的环节0和负值会让这个环节崩溃。解决办法是检查性状数据如果是测量上有真0值比如无种子要想想这个0是真实测量还是缺失值缺失了就不要填0应该填NA或者用插补。二类是提示Species present in the community matrix are missing from the trait matrix——这说明物种名没对上回到上一节说的数据体检步骤。3.3 批量计算多群落功能多样性并导出结果dbFD()的好处是一次性能算出所有样方的指数不需要循环。结果是一个列表要提取每个指数并整理成数据框。下面是我常用的提取代码# 提取功能多样性指数为数据框 fdiv_df - data.frame( Site rownames(comm), FRic fdiv_out$FRic, FEve fdiv_out$FEve, FDiv fdiv_out$FDiv, FDis fdiv_out$FDis, RaoQ fdiv_out$RaoQ ) fdiv_df这样就得到了一个行是样方、列是指数的数据框可以直接和环境的变量做后续分析比如多元回归、排序、方差分析。如果要导出成CSVwrite.csv(fdiv_df, functional_diversity_indices.csv, row.names FALSE)导出之前建议用一个检查函数抓一下数据合理性至少确认每个指数是否在合理值域内。FRic的取值范围取决于性状量纲和物种数不好做通用判断FEve和FDiv理论上是0到1之间FDis和RaoQ的取值范围不是0到1它们受性状距离尺度影响超过1很正常别看到大于1就以为算错了。3.4 关于零多度物种的处理很多人忽略的一个问题如果一个物种只在少数样方里出现在多数样方里多度为0但这个物种的性状数据仍然参与功能距离矩阵的构建。这样计算时这个物种会在那些没有它出现的样方之外提供一个隐性功能背景。FD包的处理方式是对每个样方只基于实际出现的物种计算凸包和距离所以理论上不影响各指数的计算结果。但如果一个物种在所有样方里都只有极低多度它的存在感非常弱不会显著改变群落的功能空间。如果你的数据里有大量罕见种我的建议是设置一个最小多度阈值比如将在所有样方中相对多度都低于1%的物种剔除这样能减少计算噪声但阈值的选择要写清楚并最好做敏感性分析——阈值设0.5%和2%分别跑一遍如果结论方向一致就问题不大。4. 常见问题与排查技巧实录4.1 报错速查与解决对照表我把自己在培训、答疑时遇到的高频报错和解决办法整理成一张速查表这些报错信息不一定完全原样但特征非常明显报错/警告信息典型原因解决办法Error in dbFD: Species present in the community matrix are missing from the trait matrix物种名不匹配或者性状数据行数少于群落矩阵列数检查物种名拼写、空格、大小写用intersect()提取共有物种Error in dbFD: The x matrix has too many columns compared to rows数据方向反了群落矩阵应该是行样方、列物种用t()转置矩阵Error in dbFD: some trait variables are not numeric性状矩阵里有因子或字符列用str(trait)检查每一列类型把分类变量转为数值或虚拟变量Error in dbFD: Zero or negative values in a matrix性状数据中有0或负值导致对数变换失败检查原始数据缺失值填NA而非0Warning: FEve is NA for some sites样方内物种数太少功能均匀度无法计算检查样方物种数确保每个样方至少3个物种Warning: Results may be unreliable if species number is low relative to trait number物种数相对性状数太少功能空间维数过多减少性状数量或对性状做主成分分析后取前几轴4.2 结果解读别只看显著性计算功能多样性指数很容易难的是解读。我见过不少人在得到FRic之后直接做单因素方差分析得出有差异的结论就完事了但三个细节值得再想想。第一FRic对物种数高度敏感。样方里的物种数越多凸包体积天然就越大。如果要作组间比较建议先检查组间物种丰富度是否均衡。如果一组因为干扰导致物种数大幅下降FRic下降是必然的这不能说明功能策略的压缩只是物种丢失的伴随效应。更严谨的做法是同时分析物种丰富度和FRic用物种数做协变量或者做稀疏化处理后比较。第二FEve和FDiv放一起看会比单看一个更有效果。FEve低说明功能空间里有些区域物种多度很高、有些区域很空FDiv高说明多度集中在极端物种上。这两个指数一组合你能判断出是多度集中在功能空间边缘的少数物种还是多度均匀铺展在整个功能空间对理解群落的资源利用策略很有帮助。第三qual.FDis这个输出很多人忽略了。它衡量的是用少数几个维度还原物种间原始距离的质量数值越接近1说明降维损失越小、功能空间还原越完整。理想情况下应该大于0.8或0.9如果远低于这个值说明你选的性状不能很好区分布物种需要考虑增加性状或者重新审视性状选择。4.3 我的数据清洗心得与工作流建议最后把我的一套工作流原样分享出来。前期花半小时整理数据效果远好于后期反复试错。第一步拿到原始数据先画物种累计曲线和样方物种数摸清数据集的基本盘。这个过程不需要功能多样性但能让你知道哪些样方要剔除。比如某个样方只有1个物种它算出来的FRic就毫无意义。第二步写一个数据体检函数自动检查群落矩阵和性状矩阵的物种名一致性、是否有重复、是否有NA值。我的代码如下check_fd_data - function(comm, trait) { # 检查物种名匹配 sp_comm - colnames(comm) sp_trait - rownames(trait) cat(群落矩阵物种数, length(sp_comm), \n) cat(性状矩阵物种数, length(sp_trait), \n) cat(共有物种数, length(intersect(sp_comm, sp_trait)), \n) # 检查NA if(any(is.na(trait))) cat(警告性状矩阵存在NA值\n) if(any(is.na(comm))) cat(警告群落矩阵存在NA值\n) }第三步多次运行dbFD()前先在文档里写下你对结果的先验预期。比如我预期干扰度高的样方应该有更低的FRic、更高的FDiv那么跑完结果拿这个预期一对照如果完全相反不是你的预期错了多半是数据格式或者物种筛选出了问题。这个预期验证策略帮我抓出过很多肉眼难发现的错误。我一再强调数据格式整理的重要性是因为这个环节真的决定成败。很多发表的结果后来被发现不可复现追溯原始分析脚本往往是数据规整环节的一个bug导致的。养成先体检后计算先预期后解读的习惯能帮你减少大量返工时间。工具本身不复杂复杂的是你有没有一套隐形的工作规范。
阅读完成 · 觉得有帮助？
