☰
R语言生态位宽度计算与可视化完整指南
2026/10/3 2:33:14 网站建设 项目流程

做生态数据分析的人,十有八九都会卡在生态位宽度这个指标上。倒不是计算本身有多难,而是市面上的教程要么只讲公式不讲实现,要么给了代码却不说数据该怎么整理、结果该怎么解读、图该怎么画才拿得出手。我这两年用R语言做了不少物种资源利用方面的分析,从最基础的Levins指数到可视化出图,踩过的坑不少,也摸索出一套完整可复用的流程。这篇就把整个链路拆开讲清楚,从数据准备到算法原理,从R代码实现到发表级图表绘制,一次性讲透。

1. 生态位宽度到底是什么

生态位宽度(niche breadth)衡量的是一个物种对资源利用的范围和多样程度。说得直白点,如果一个物种只在一种资源上活动,那它的生态位就很窄;如果它能在多种资源之间灵活切换,生态位就宽。这个指标在物种共存机制、群落构建、入侵生态学、食性分析这些研究方向里都是高频出现的核心参数。

计算生态位宽度的核心输入是一张“物种-资源利用矩阵”,行是物种,列是资源状态,单元格里的数值可以是多度、频次、生物量占比,也可以是取食记录数、环境样本检出数。矩阵准备好之后,常用算法有三类:Levins、Shannon-Wiener和Hurlbert。这三种算法的思路差异很大,选哪个取决于你的数据形态和研究问题。

Levins指数是最经典的做法,它直接衡量资源利用的均匀程度。公式是 B = 1 / Σ(pᵢ²),其中pᵢ是第i种资源利用比例,也就是该物种在资源i上的利用量占它总利用量的比例。这个指数的直观含义是:一个物种的利用比例越集中在少数资源上,Σ(pᵢ²)就越大,B就越小;利用越分散均匀,B越大。它还有一个标准化版本 B_ST = (B - 1) / (n - 1),n是资源总数,标准化后范围落在0-1之间,便于跨研究比较。

Shannon-Wiener指数则是从信息论角度切入,公式是 H' = -Σ(pᵢ × ln pᵢ)。它同时对资源利用的丰富度和均匀度敏感,数值越大代表生态位越宽。这个指标在生态学里用得非常多,但要注意它受样本量和资源分类粒度的影响较大,横向对比时要谨慎。

Hurlbert指数解决了一个实际问题:资源本身的稀缺或丰富程度不一样,如果一个物种利用了某种常见资源,和利用了某种稀有资源,对生态位宽度的贡献理应不同。Hurlbert在公式里引入了资源丰度权重,更贴近真实生态情景。

为了让你直观感受三种算法的差异,我跑了一组模拟数据:5个物种、6种资源的利用矩阵,计算得到如下结果。

物种Levins BLevins B_STShannon H'Hurlbert B_H
物种A1.910.180.861.42
物种B3.750.551.452.80
物种C1.260.050.420.95
物种D4.820.761.683.66
物种E2.630.331.121.98

物种D在所有资源上利用都相对均衡,所以四种指数都给出最高值;物种C几乎只在一种资源上活动,所有指数都最低。但注意排序上并不完全一致,这就是因为不同算法对“宽度”的定义侧重不同。实际发表文章时,通常选一种主指标,再辅以另一种做敏感性验证,而不是把所有指数都堆上去。

2. 数据准备与预处理

很多人在R里跑生态位宽度计算,数据格式就卡住了。常见的数据来源是样方调查表,比如每一行是一个样方、每一列是一个物种的多度,或者反过来——每一行是一个物种、每一列是一种资源类型。做生态位宽度分析,必须先把数据整理成“物种为行、资源为列”的矩阵格式。

举个例子,假设你调查了5个物种在6种资源上的利用频次,数据框至少要长这样:

speciesres1res2res3res4res5res6
sp11205380
sp20152094
sp36301000
sp4127456
sp5300001

注意几点:行名必须是物种名,不能是数字序号,否则后面计算时标签会丢失;列名是资源状态,建议用带单位或明确含义的名称,方便出图时直接作为坐标轴标签;0值表示未利用,这是正常情况,不要删掉。我自己习惯把所有数值列保持为numeric类型,如果从Excel复制进来变成了character,直接用mutate(across(where(is.character), as.numeric))批量转换。

还有一种常见情况:数据是长格式的,每一行是一个观测记录,包含物种名和资源类型两列。这种格式不能直接用于生态位宽度计算,需要先转化为宽格式。用tidyr包的pivot_wider即可。比如原始数据有三列:species、resource、count,转换代码是:

library(tidyr) wide_data <- raw_data %>% pivot_wider(names_from = resource, values_from = count, values_fill = 0) %>% column_to_rownames("species")

这里values_fill = 0很重要,它的作用是把没有记录的组合补0,避免生成NA导致后续计算报错。长转宽之后,再确认一下矩阵维度,行数等于物种数,列数等于资源数,就可以正式进入分析了。

接下来是标准化问题。如果你的数据是不同采样强度下汇总的,比如有的样方调查面积大、有的调查天数多,那最好先做标准化。最简单的标准化是把每个物种的利用量除以它的总和,得到资源利用比例,也就是pᵢ矩阵。但要注意,如果某些物种的总利用量极小,标准化后比例会显得很分散,解释时要小心。

还有一种情况是资源本身的丰度存在巨大差异。比如食性分析中,有些猎物种类在环境里本身就很丰富,有些很稀有。这时候用Hurlbert指数比Levins更合适,因为Hurlbert能把这些背景差异纳入考量。如果研究设计里没有环境资源丰度数据,那就老老实实用Levins或Shannon,并在方法部分说明理由。

预处理阶段我还会顺手做几件小事:检查有没有全为0的行(说明该物种没有利用记录,应删除或调查确认),检查有没有全为0的列(说明该资源没有被任何物种利用,删除后能避免分母虚高),以及看有没有极端离群值。生态位宽度对极端值敏感,如果某个单元格的值是其他值的几十倍,建议先确认不是录入错误。

3. R语言计算生态位宽度的三种实现方式

3.1 用spaa包一键计算

生态学分析圈里,spaa这个包就是为这类问题设计的。它专门处理物种生态位和生态位重叠计算,函数接口很简洁。安装和加载方式:

install.packages("spaa") library(spaa)

核心函数是niche.width(),它默认计算Levins和Shannon两个指数。直接把宽格式矩阵喂进去就行:

niche_result <- niche.width(wide_data, method = "levins") niche_result <- niche.width(wide_data, method = "shannon")

但要注意,spaa包要求输入必须是矩阵格式,数据框会报错。杠精一点的读者可能已经发现了:明明as.matrix转一下就能解决。没错,我在调包之前都会加一个as.matrix(),并且把行名列名都确认清楚。还有一个坑是method参数的取值——"levins"和"shannon"是它认识的两个写法,不区分大小写的事我没实测过,但建议严格用小写。

如果你需要Hurlbert指数,spaa包默认不带这个功能,需要自己写函数。下面给出我一直在用的版本,参数里多了一个resource.freq,用来传入各资源在环境中的丰度比例向量:

hurlbert_width <- function(mat, resource.freq = NULL) { if (is.null(resource.freq)) { resource.freq <- rep(1, ncol(mat)) } else { resource.freq <- resource.freq / sum(resource.freq) } p <- mat / rowSums(mat) # 计算每个物种的Hurlbert生态位宽度 B_H <- apply(p, 1, function(x) { 1 / sum((x^2) / resource.freq) - 1 }) return(B_H) }

这段代码的逻辑是:把利用比例pᵢ除以资源相对丰度qᵢ,再平方求和取倒数,最后减1。减1的目的是让最小可能值变为0,方便解释。实际使用中记得传入的resource.freq向量长度必须等于矩阵列数,否则R会静默按循环规则补齐,结果完全错掉。

3.2 手写公式加深理解

如果你不想只做调包侠,手写一遍公式能让你对算法理解深很多。生态位宽度计算本质上是矩阵运算和向量运算的组合,代码量并不大。

Levins生态位宽度的手写实现:

levins_b <- function(mat) { p <- mat / rowSums(mat) B <- 1 / rowSums(p^2) B_ST <- (B - 1) / (ncol(mat) - 1) return(data.frame(B = B, B_ST = B_ST)) }

Shannon-Wiener指数的手写实现:

shannon_h <- function(mat) { p <- mat / rowSums(mat) H <- -rowSums(p * log(p), na.rm = TRUE) return(H) }

写的时候有个小细节容易出错:p * log(p)在p为0时会得到NaN,因为0乘以负无穷。所以要么先对p做过滤,要么在rowSums里加na.rm = TRUE。我习惯用后者,简单粗暴。

手写的好处是方便自定义。比如有的生态学文献里会对Shannon指数做e的幂次变换,得到有效物种数形态的宽度值;有的会要求基数为2而不是自然对数。这些都可以在公式层面灵活修改,而不需要去改包源码。

3.3 批量计算多组数据

实际分析中经常要对多个群落或者多个年份的数据分别计算生态位宽度。这时候与其手动循环,不如把数据整理成嵌套列表或者带分组信息的长表,然后用purrr包批量处理。

假设你有一个community列表示不同样地或年份,数据已经按物种-资源矩阵分组存放,可以用嵌套数据框的方式批量跑:

library(purrr) library(dplyr) result <- raw_data %>% group_by(community) %>% nest() %>% mutate(niche = map(data, function(df) { mat <- df %>% column_to_rownames("species") %>% as.matrix() niche.width(mat, method = "levins") })) %>% unnest(niche)

这个流程的核心是nest()把每个分组的矩阵变成一个独立的数据框,再用map逐一计算,最后unnest展开结果。输出就是一个包含群落名、物种名和生态位宽度值的整洁数据框,后续无论是做统计检验还是画图都非常方便。

多说一句:nest()和unnest()这对函数是tidyr包在1.0版本之后的标志性功能,早期版本没有。如果你装了老版本R,记得先install.packages("tidyr")升级。

4. 生态位宽度的可视化展示

计算完生态位宽度只是第一步,能画出清晰的图才算真的做完分析。可视化不光是给论文配图用的,我自己在做探索性分析时也会画一堆图来快速发现模式。下面按从简单到复杂的顺序,讲几个亲测好用的方案。

4.1 基础条形图:快速对比物种间宽度

条形图是生态位宽度展示最直接的方案,适合物种数量在10-30个范围内的展示。用ggplot2画,几行代码就能出图:

library(ggplot2) plot_data <- result %>% rownames_to_column("species") ggplot(plot_data, aes(x = reorder(species, B), y = B)) + geom_col(fill = "#4E79A7", width = 0.7) + coord_flip() + labs(x = "物种", y = "Levins生态位宽度 (B)") + theme_minimal(base_size = 14) + theme(panel.grid.major.y = element_blank())

这里有两个技巧值得留意。一是用reorder(species, B)按宽度值排序,画出来的图呈现从左到右递增或递减的阶梯效果,比随机排序好看得多。二是用coord_flip()把条形变成横向——物种名一般比较长,横向排布能完整显示标签,避免文字互相挤压。

如果你要在一张图里同时展示多个指数(比如Levins和Shannon),建议用facet_wrap分面,而不是硬塞进同一坐标系,因为不同指数的数值范围差异很大,放在一个面板里会有视觉误导。

4.2 雷达图:展示资源利用谱

条形图只能看宽度值,看不到物种到底用了哪些资源。雷达图(蜘蛛网图)可以把每个物种的资源利用谱完整展示出来,特别适合物种数量少(3-8种)、资源维度明确(比如食物类型、栖息地类型)的场景。R里画雷达图我用的是fmsb包:

install.packages("fmsb") library(fmsb) # 准备数据:行为资源,列为物种 radar_data <- t(as.data.frame(p)) radar_data <- as.data.frame(radar_data) # fmsb要求前两行是最大值和最小值 radar_data <- rbind(rep(max(radar_data), ncol(radar_data)), rep(0, ncol(radar_data)), radar_data) radarchart(radar_data, axistype = 1, pcol = c("#E64B35", "#4DBBD5", "#00A087", "#3C5488", "#F39B7F"), pfcol = scales::alpha(c("#E64B35", "#4DBBD5", "#00A087", "#3C5488", "#F39B7F"), 0.25), plwd = 2, cglcol = "grey60", cglty = 1, axislabcol = "grey30", vlcex = 0.8) legend(x = 1.2, y = 1.1, legend = rownames(radar_data)[3:nrow(radar_data)], bty = "n", pch = 20, col = c("#E64B35", "#4DBBD5", "#00A087", "#3C5488", "#F39B7F"), text.col = "black", cex = 0.9, pt.cex = 1.5)

雷达图有一点反直觉:面积大不代表生态位宽,因为某些资源维度跨度大(比如从浅水到深水)会拉伸图形。所以看雷达图时重点看形状,而不是面积。如果所有物种的雷达图叠在一起太乱,分面画才是正解,fmsb包不支持分面,配合gridExtra包手动排版可以解决。

4.3 热图:同时展示多个物种和多种资源

热图能把资源利用矩阵本身的模式可视化,本质上和生态位宽度相辅相成。它展示的是pᵢ矩阵,也就是每个物种在每种资源上的利用比例,颜色深浅代表比例高低。这样能直观看出哪些物种是“专性利用”(图上一片深色集中在一个角落),哪些是“广谱利用”(颜色分布均匀)。

ggplot2+geom_tile()是最清晰的热图方案。配合scale_fill_gradient2()做发散色标,0用白色表示,高值用暖色,效果很好:

p_long <- p %>% as.data.frame() %>% rownames_to_column("species") %>% pivot_longer(-species, names_to = "resource", values_to = "prop") ggplot(p_long, aes(x = resource, y = species, fill = prop)) + geom_tile(color = "white", linewidth = 0.5) + scale_fill_gradient2(low = "white", mid = "#FFD36E", high = "#C70039", midpoint = 0.3, limits = c(0, 1)) + labs(x = "资源类型", y = "物种", fill = "利用比例") + theme_minimal(base_size = 13) + theme(axis.text.x = element_text(angle = 45, hjust = 1))

热图的排序很关键。如果按字母序排列物种和资源,画出来的图就是一片噪声。我通常先用层次聚类对行列重排,让相似的物种和资源靠在一起,图案才有规律可寻。pheatmap包一步到位支持行列聚类,推荐试试——在生态位分析里算是潜力被低估的工具。

4.4 排序图:把生态位关系放到群落语境里

更进阶的做法是把生态位宽度放回到群落排序(ordination)的语境中来解读。常用方案是先做CA或CCA分析,再把每个物种的生态位宽度映射到排序图的点大小或颜色上。这样既能看物种在资源轴上的位置,又能看它利用范围的大小。

vegan包是排序分析的标准工具。示例:

library(vegan) # 对利用矩阵做对应分析 ca_model <- cca(wide_data) # 提取样方和物种坐标 site_scores <- scores(ca_model, display = "sites") species_scores <- scores(ca_model, display = "species") # 用生态位宽度作为点的大小 species_df <- as.data.frame(species_scores) species_df$B <- niche_result$B ggplot(species_df, aes(x = CA1, y = CA2, size = B)) + geom_point(color = "#2E86AB", alpha = 0.8) + geom_text(aes(label = rownames(species_df)), vjust = -0.8, size = 3) + scale_size_continuous(range = c(2, 10)) + labs(x = "CA1", y = "CA2", size = "生态位宽度") + theme_minimal(base_size = 14)

这种图在论文里特别讨喜,因为一箭双雕:既展示了物种在资源空间中的分化,又突出了宽度差异。如果你的研究重点是群落构建机制,这张图比单独的柱状图有信息量得多。

5. 常见问题与排查技巧

5.1 为什么计算结果和文献对不上

这是被问得最多的问题。DNA条形码研究里,同样一个物种,不同文献报道的生态位宽度差异巨大,误差来源往往不是算法错了,而是定义不同。有人用的是标准化Levins(B_ST),有人用的是原始B值;有人对数据先做了平方根变换,有人是原始多度。另外,资源分类粒度不同也会导致生态位宽度明显不同——同样是食性分析,猎物按“目”分和按“种”分,算出来的宽度当然不可比。建议分析前先把资源分类方案写清楚,并在方法部分明确列出用的是哪个指数、哪个公式版本。

5.2 出现0值导致计算结果异常

如果某个物种在矩阵里全是0,rowSums就是0,除以0得到NaN。更隐蔽的问题是某个物种只利用了一种资源,这时候p矩阵里有一个1和一堆0,Levins B正好等于1,B_ST等于0,看起来正常,但如果某种资源的列和是0,会让某些物种的pᵢ出现NaN。处理方式是在计算前删除全零行和全零列。

5.3 可视化时中文标签乱码

这个坑在Windows系统上特别常见。R的默认绘图设备对中文支持不友好,ggplot2画图时中文会显示为方格。解决方案有两个:一是安装中文字体并设置theme(text = element_text(family = "STHeiti"))之类的参数,二是干脆把图内标签统一用英文,等投稿时再在图片处理软件里加中文注释。我个人的习惯是出图时用英文,在论文正文用中文描述,省去很多字体兼容问题。

5.4 样本量不均衡导致比较失真

生态位宽度受资源取样充分度影响很大。一个物种调查了200条取食记录,另一个只有15条,前者的生态位宽度天然会更大。这个问题没有完美的统计学修正,但有一个靠谱的应对策略:用稀松曲线(rarefaction curve)检查采样充分性,或者对数据进行多次重抽样后计算宽度的置信区间。vegan包里的rarecurve()函数可以用于检查样本量是否足够。

问题症状解决方案
数据格式错误niche.width()报错用as.matrix()转换,确保行名为物种名
全零行/列结果出现NaN或Inf删除全零行/列,再重新计算
中文乱码出图中文显示为方格设置中文字体或改用英文标签
三种指数结果排序不一致物种间宽度排名变化以研究问题为核心选主指标,其他做敏感性分析
样本量差异大宽度值虚高或虚低做稀松曲线检查采样充分性,或重抽样计算置信区间

6. 分析报告与成果输出建议

生态位宽度计算和可视化做完之后,还有一个容易被忽视的环节:如何把分析流程和结果完整输出。我见过不少人在R Studio里运行完代码,窗口一关,三个月后要写论文时完全想不起来当时做了什么处理。这里给两个实用建议。

一是用R Markdown记录整个分析流程。R Markdown能让代码、结果、图表和文字说明写在同一个文档里,输出成HTML或PDF就是一份完整的分析报告。生态位宽度分析这种流程相对固定的任务,非常适合做成模板,下次换一批数据直接替换矩阵内容就能复用。

我自己的模板结构是这样的:

--- title: "生态位宽度分析报告" output: html_document --- ## 数据说明 - 数据来源:XXX调查 - 资源分类方案:XXX - 采样时间:XXX ## 数据预处理 ```{r, echo=TRUE} # 读取和整理数据

生态位宽度计算

# 计算Levins和Shannon指数

可视化

# 条形图和热图

结果解读

  • 生态位最宽的物种:XXX
  • 生态位最窄的物种:XXX
二是导出所有数值结果时,最好保留足够多的小数位数,并同时导出原始p矩阵和最终宽度值。后续如果审稿人要求换一种指数重算,或者要求补充Bootstrap置信区间,你不需要重新跑整个流程,直接修改一小段计算代码就行。 这里再分享一个我常用的经验:生态位宽度的Bootstrap置信区间。这个方法能帮你评估宽度估计的稳定性,特别是样本量不足的时候非常有用。做法是对每个物种的利用记录做有放回重抽样(次数=原样本量),重复1000次,每次计算宽度值,最后取2.5%和97.5%分位数。代码实现基于`boot`包或者自己写`for`循环都可以,我自己更常用`purrr::map_dbl`配合`sample`来做。 ```r boot_niche_width <- function(vec, n_boot = 1000) { boot_vals <- map_dbl(1:n_boot, function(i) { boot_vec <- sample(vec, replace = TRUE) p <- boot_vec / sum(boot_vec) 1 / sum(p^2) }) quantile(boot_vals, c(0.025, 0.975)) } # 示例:对物种A的利用频次做Bootstrap vec_A <- as.numeric(wide_data["sp1", ]) boot_niche_width(vec_A)

不过要提醒一点:Bootstrap适用于利用记录是“独立取样”的场景,如果数据本身是空间自相关的,这个方法会产生偏乐观(过窄)的置信区间。这种情况下更严谨的做法是分层Bootstrap或空间块Bootstrap,但生态学论文里能做到第一种就已经算考究了。

关于生态位宽度计算,我实际摸爬滚打这么久,最大的感受就是:计算本身不是难点,难的是对指标含义的理解和对数据质量的把控。很多初学者拿到数据就急着跑代码,等到画图的时候才发现前期的数据整理有问题,返工成本极高。我现在的习惯是:先花40%的时间把数据格式、标准化方案、算法选择理清楚,再花20%时间写代码,剩下40%时间都在做可视化调优和结果解释。这个过程急不得。如果你正准备开始做类似的分析,建议从一张干净的物种-资源矩阵开始,先跑通流程,再逐步增加复杂度,比如加入环境因子、空间信息或者多个群落的比较。生态位宽度是一个非常好的切入点,弄懂它之后,生态位重叠、生态位分化这些概念都会一通百通。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询