R语言生态位宽度计算实战：从Levins到Shannon的完整指南

前阵子帮一位做农学研究的师弟处理群落数据，他拿着一张物种-生境记录矩阵来找我，说文献里同一种天敌昆虫的生态位宽度，不同文章写出来的数值差得离谱，自己在R语言里算的Levins指数跟手算也对不上。这个问题看着简单，其实牵出一连串容易翻车的细节：指标选哪个、数据要不要归一化、标准化公式里那个“-1”是干什么的、画出来的图怎么才能有信息量而不是花架子。我在这篇文章里把完整的计算流程、可视化方案和踩坑记录都梳理一遍，适合生态学、农学、保护生物学方向的研究生，以及任何需要用R语言做群落数据分析的人。

1. 生态位宽度到底是什么，为什么我总是绕不开它

1.1 一个听起来简单、做起来容易翻车的概念

生态位宽度（niche breadth）衡量的是一个物种对资源利用的多样化程度。说人话就是：这个物种到底是“什么都吃”的广谱选手，还是“只认准一种”的专性选手。广义种像小区门口什么菜都卖的饭馆，狭食种像只做一道招牌菜的苍蝇馆子，各有各的生存策略。

这个指标在群落生态学里几乎绕不开，因为它直接关系到三个核心问题：种间竞争是否激烈、群落结构是否稳定、环境变化来临时哪些物种更容易受冲击。具体到天敌昆虫，宽生态位的捕食者往往能在多种农田生境间移动，对害虫种群的控制也更持久；而窄生态位的物种虽然敏感脆弱，但通常是特定生境的指示生物。

做这类分析时，R语言是我最顺手的工具，因为从数据清洗、指标计算到可视化，一套tidyverse流程就能全部跑通。不过先用哪个指标、怎么标准化、数据从长表转成宽表时有没有坑，这些细节如果不提前想清楚，后面算出来的数字很容易让人误读。

1.2 资源状态怎么划分，决定后续所有计算的边界

“资源状态”这个词听起来很学术，但它可以非常接地气：食性分析中的猎物种类、栖息地选择中的生境类型、昼夜活动节律中的时间段、传粉网络中的访花植物类别，都可以作为资源状态。

资源状态的划分有两条硬性要求：第一是互斥，一个记录只能落在一个状态里，比如“稻田”和“路边”如果在地理上有重叠，那就要重新定义类别；第二是可识别，每个状态在野外或实验设计里必须是能稳定区分的。

经常有人在这里栽跟头——把连续环境变量（比如土壤含水量0%-100%）直接切成等距区间，切几段？每段代表什么生态学含义？如果切得太粗，宽度会被系统性高估；切得太细，很多物种的记录数不够，p_i全是零碎小数，计算不稳定。我的建议是先做频次直方图，看记录在资源轴上的实际分布再决定切分边界，而不是机械地等分。

2. 三种主流指标的计算逻辑与选型指南

2.1 Levins宽度：均匀度视角

Levins（1968）提出的公式是应用最广泛的：

B = 1 / Σ(p_i²)

其中p_i是该物种对第i个资源状态的利用比例。这个公式的直觉很直接：如果只利用一种资源，p只有一个1，其他全是0，B=1，宽度最小；如果完全均匀地利用n种资源，每个p_i都等于1/n，B=n，宽度最大。

但直接用B有个问题：资源状态数n不同的时候，B的数值上限不一样。同样是B=3，在一个5状态系统里是中等宽度，在一个10状态系统里就偏窄了。所以实际工作中更常用标准化形式：

Ba = (B - 1) / (n - 1)

标准化后取值在0到1之间，方便跨数据集比较。这个“-1”和“n-1”不是数学炫技，而是把B的起点压到0、最大值压到1。

2.2 Shannon宽度：信息熵视角

Shannon宽度直接借用信息论里的熵：

H = -Σ(p_i × ln(p_i))

它衡量的是“猜这个物种下一次会出现在哪个资源状态时的不确定性”。利用比例越均匀，不确定性越大，H越高。取值范围在0到ln(n)之间，标准化方式为H / ln(n)。

Shannon指标对稀有资源状态的“贡献”比Levins更敏感，因为对数项会让那些利用比例很低的末端资源也不被直接忽略。如果你关心的是物种对边缘性资源的利用潜力，Shannon会比Levins更合适。而且它的计算方式跟α多样性里的Shannon指数一模一样，很多人会把这一套逻辑顺延过来，学习成本很低。

2.3 Smith与Hurlbert：把资源可得性放进来

Levins和Shannon都默认所有资源状态在环境中出现机会相同，这在野外往往不成立。如果农田里麦田面积占比60%、果园只占5%，那一个物种在麦田里记录多，可能仅仅是因为麦田更容易被碰到，而不是它真的偏好麦田。

Smith（1982）提出：

FT = Σ√(p_i × a_i)

其中a_i是第i种资源在环境中的可得比例。这个指标衡量的是物种实际利用模式与资源可得模式的匹配程度，取值0到1，越接近1表示利用比例与可得比例越一致。

Hurlbert（1978）的思路类似：

B = 1 / Σ(p_i² / a_i)

它相当于给每个p_i都除以对应资源的可得性，把“使用多”修正为“相对于可得性而言使用多”。这两个指标都需要额外收集环境资源比例数据，不是所有研究都具备条件，但金标准意义上它们比Levins和Shannon更贴近生态学现实。

2.4 指标选型对照表

选型的核心原则：手上只有利用频次数据时，优先Levins标准化值或Shannon标准化值；如果研究设计里已经测了资源可得性，就不要再回避Smith或Hurlbert，否则审稿人大概率会问。

3. 从原始记录到干净矩阵：数据准备的核心细节

3.1 长表转宽表：pivot_wider的大坑

野外调查数据最常见的存储格式是长表：每一行是一个样本记录，列包括调查点、物种名、资源状态、记录数。计算生态位宽度时需要把它转成宽表矩阵，行是物种，列是资源状态，单元格是记录数或比例。

tidyverse里的pivot_wider是标准做法：

library(tidyverse) raw_data <- read_csv("field_records.csv") mat <- raw_data |> group_by(species, resource_state) |> summarise(count = sum(count), .groups = "drop") |> pivot_wider(names_from = resource_state, values_from = count, values_fill = 0) |> column_to_rownames("species") |> as.matrix()

这里最容易犯的错误是忘记values_fill = 0。野外记录里没有出现过的组合通常是空行，pivot_wider默认会填成NA，如果不补成0，后面sum计算会一并把NA卷进去，生态位宽度直接变成NA。

3.2 0、NA和缺失值

生态学数据里的0和NA含义完全不同。0代表“调查了但没记录到”，是有信息量的数据点；NA代表“没有调查”或“数据丢失”，是不能参与计算的。

所以宽表矩阵里填0是安全的，但要注意区分真正的缺测。如果把一个根本没调查的生境类型填成NA，Levins计算时会把整行当作缺失处理，结果全组物种都报错。遇到这种情况，要么删除该资源状态，要么用多重插补或半定量估计补齐，绝对不能留NA进公式。

另外我建议在矩阵生成后先做一次全面检查：

summary(mat) any(is.na(mat))

如果检查出NA，优先回溯原始记录确认是“没调查”还是“没记录到”，前者走删除列方案，后者填0。

3.3 抽样强度不一致怎么办

这是生态位宽度计算里最隐蔽的系统性误差。Levins和Shannon的公式本身对总记录量做了归一化，所以一个物种记录100条和记录300条，计算出来的宽度数值不在同一个统计功效水平上。

比如A物种只被调查到15条记录，恰好集中在两个生境里，算出来宽度很窄；B物种被系统调查了500条，覆盖五个生境，算出来宽度很宽。这个对比是站不住的，前者很可能只是采样不够。

解决办法有三个层次：最理想是原始调查时就做均匀抽样设计；已经拿到数据的，可以做稀释（rarefaction），用vegan包的rrarefy把各物种记录量抽到同一水平再算；实在不能抽稀的，至少要严格控制数据量差异，并报告各物种的总记录数，让读者自己判断可信度。很多“花里胡哨”的生态位宽度论文，问题就出在这一步。

4. R语言实现：手写函数与spaa包双方案对照

4.1 手写Levins与Shannon计算函数

自己写函数的好处是逻辑透明，出了问题一眼就能定位。下面这段代码可以直接贴进RStudio运行：

# Levins 生态位宽度 levins <- function(x) { p <- x / sum(x) B <- 1 / sum(p^2) n <- length(x) Ba <- (B - 1) / (n - 1) c(B = B, Ba = Ba) } # Shannon 生态位宽度 shannon_width <- function(x) { p <- x / sum(x) p <- p[p > 0] # 去掉0，避免log(0) H <- -sum(p * log(p)) n <- length(x) Hstd <- H / log(n) c(H = H, Hstd = Hstd) } # 逐行计算 apply(mat, 1, levins) apply(mat, 1, shannon_width)

apply(mat, 1, levins)会返回一个两行矩阵，第一行是B，第二行是Ba，列名对应物种名。Shannon同理。

这里有个细节：shannon_width里p <- p[p > 0]不能省略。如果一个物种在某资源状态上没有记录，p=0，log(0)直接返回-Inf，整行结果全是NaN。

4.2 更高阶：用spaa包一个函数跑完三套指标

如果你不想手写，生态学常用的spaa包里有现成函数niche.width，它能同时计算Levins、Shannon和Smith三类指标：

library(spaa) # mat：行是物种，列是资源状态 res_levins <- niche.width(mat, method = "levins") res_shannon <- niche.width(mat, method = "shannon") # Smith指标需要额外指定资源可得比例向量 a a <- c(0.20, 0.25, 0.20, 0.15, 0.20) # 假设五类生境的面积占比 res_smith <- niche.width(mat, method = "smith", A = a)

用这个包前我会建议先跑一遍str(res_levins)看输出结构。不同版本包里返回对象可能是矩阵也可能是数据框，列名有时是B和Ba，有时是Levins和Levins.std，直接print容易看漏。记住这个原则：任何R包返回的复杂对象，先str()，再取数。

4.3 双方案对拍与结果规整

手写函数和R包的结果必须保持一致，这一步叫“对拍”。我在第一次用spaa时发现手写版和包版本差了小数点后几位，后来发现是Smith指标里sqrt(0)的浮点数精度问题，不是算法错误，对生态学结论没有影响。

不管用哪种方案，最后都要整理成一个规整的数据框，方便join和画图：

result_df <- data.frame( species = rownames(mat), Ba = res_levins[, "Ba"], Hstd = res_shannon[, "Hstd"] ) result_df |> arrange(desc(Ba))

5. 可视化不是画图而已：四类图表的表达逻辑

5.1 条形图：物种间的宽度排序

计算完一堆数字之后，第一张图应该是物种间宽度的排序条形图。关键点在于“排序”而不是随便画：宽度值只有横向对比才有意义，不排序的条形图信息量直接减半。

library(ggplot2) result_df |> mutate(species = fct_reorder(species, Ba)) |> ggplot(aes(x = species, y = Ba)) + geom_col(fill = "#4C72B0", width = 0.6) + coord_flip() + labs(x = NULL, y = "标准化Levins生态位宽度") + theme_minimal(base_size = 13)

从这张图能很快看出谁是大范围活动者、谁高度特化。给论文配图时，我习惯再叠加一个Shannon标准化结果的副面板，两列对比，因为Levins和Shannon排序结果有时会不同，这种差异本身就很有生态学故事可以讲。

5.2 资源利用曲线：形状比数值更会说故事

生态位宽度只是一个综合数值，它掩盖了资源利用曲线的具体形状。两个物种可能计算出的宽度完全一样，但一个偏嗜两三种资源，另一个对所有资源雨露均沾，生态学含义截然不同。

把利用比例画成折线图就能看到这些模式：

prop_df <- mat |> as.data.frame() |> rownames_to_column("species") |> pivot_longer(-species, names_to = "habitat", values_to = "count") |> group_by(species) |> mutate(prop = count / sum(count)) |> ungroup() ggplot(prop_df, aes(x = habitat, y = prop, color = species, group = species)) + geom_line(linewidth = 1) + geom_point(size = 2) + theme_minimal(base_size = 13) + labs(x = "生境类型", y = "利用比例", color = "物种")

曲线平坦的是广布型，曲线陡峭的是偏好型，曲线有几个峰的可能是资源分割型。我在实际分析中更倾向于看这张图，而不是只盯宽度数值。

5.3 热力图：完整利用格局一览

当物种数量超过10个时，折线图会变成一团乱麻，这时候热力图是最清晰的替代方案。行是物种，列是资源状态，颜色深浅代表利用比例：

ggplot(prop_df, aes(x = habitat, y = species, fill = prop)) + geom_tile(color = "white", linewidth = 0.5) + scale_fill_viridis_c(option = "C", name = "利用比例") + theme_minimal(base_size = 13) + labs(x = NULL, y = NULL) + theme(axis.text.x = element_text(angle = 45, hjust = 1))

热力图能同时看到三件事：哪些资源被普遍利用（整列颜色都很深）、哪些物种是专一性利用者（单格特别深）、哪些物种的资源谱很宽（整行颜色均匀分布）。它是最适合放进组会PPT里的图。

5.4 进阶：结合排序轴展示生态位分化

如果你的研究还采集了环境变量，可以进一步用vegan包的排序分析把物种放在生态位空间里展示。比如做RDA或NMDS，把物种点和资源变量点画在同一个双序图里：

library(vegan) # env是生境的环境变量矩阵，格式为行=调查生境，列=变量 # mat是物种-生境频次矩阵 rda_result <- rda(mat ~ ., data = env) plot(rda_result, display = c("species", "bp"))

有了排序图，物种间的宽度差异和资源利用偏向就变成了空间距离和方向，适合回答“哪些物种占了生态空间的哪个角落”这类问题。这一节属于加分项，能用上就说明你的数据质量已经过关了。

6. 完整案例：农田天敌群落的生态位宽度分析全流程

6.1 数据背景与探索目标

用一份模拟的农田天敌调查数据跑一遍全流程。数据是7种天敌昆虫在5类生境中的调查记录数：

mat <- matrix(c( 30, 45, 20, 15, 10, 5, 25, 40, 20, 10, 10, 10, 10, 10, 60, 2, 8, 5, 50, 35, 60, 10, 5, 2, 3, 25, 20, 15, 20, 20, 5, 5, 5, 5, 5 ), nrow = 7, byrow = TRUE) rownames(mat) <- c("瓢虫", "草蛉", "寄生蜂", "食蚜蝇", "步甲", "蜘蛛", "螳螂") colnames(mat) <- c("稻田", "麦田", "玉米田", "菜地", "果园")

矩阵的行是物种、列是生境类型，单元格是调查到的个体数量或记录频次。目标：计算每个物种的标准化生态位宽度，判断哪些天敌是景观尺度的广布种，哪些是局部生境的专性种。

6.2 指标结果解读

用前面的函数跑一遍后会得到类似这样的结果（不同版本浮点数略有差异）：

螳螂在五类生境中完全均匀分布，标准化宽度达到最大值1，这是典型的广布机会主义种。蜘蛛也非常接近均匀，说明它对农田景观异质性的容忍度很高。步甲则几乎锁定稻田，宽度只有0.18，专性极强，很可能和稻田湿润微环境或特定猎物有关。

寄生蜂的原始记录里60%集中在果园，宽度不高，说明它更偏好果园资源链，这可能与果园里蚜虫或介壳虫等寄主密度有关。这些解读如果不结合资源利用曲线，单看表格数字很难发现“偏好在果园”这一层信息，所以在论文里我通常会把宽度表和图2的曲线一起放。

6.3 多图联动的展示技巧

如果你的报告或论文有一个主图配额，我建议用组合图：左侧放标准化宽度条形图，右侧放资源利用曲线，下面放热力图。这样读者能从“谁宽谁窄”到“在哪个资源上宽”再到“整体格局长什么样”，逐层拆解，信息量比单张图大得多。

R里可以用patchwork包快速拼图：

library(patchwork) p_bar + p_curve + p_heat + plot_layout(ncol = 1)

当然，如果目标期刊要求矢量图，记得用ggsave把每张子图单独导出成PDF，再用AI或Inkscape进行最后排版。拼图用patchwork只是为了汇报和组会时快速预览。

7. 实战中的六个大坑和我的处理方式

7.1 坑一：把未归一化的计数直接丢进公式

有人直接用Σ(x_i²)计算，忘了除以总数得到比例。这样算出来的宽度会随总记录量剧烈变化——同一个物种调查100条算出来是45，调查1000条算出来是4500，完全失去可比性。记住：公式里的p必须是比例，所有进公式的计数都要先除行和。

7.2 坑二：Shannon指数撞上log(0)

shannon计算时如果保留了p=0的状态，log(0)会产生-Inf，连带后面的求和变成NaN。我在第一版脚本里就栽过这个跟头。解决方法是p <- p[p > 0]，只在非零比例上求和。要注意这样处理后，n在标准化公式里仍然用资源状态总数，而不是过滤后的非零个数，否则标准化结果会被错误抬高。

7.3 坑三：n=1时标准化公式直接爆炸

如果你的研究里某个资源状态轴只有1个类别（比如只调查了一种生境），(B-1)/(n-1)会变成除以0，输出Inf或NaN。这种情况在文献里也不少见——资源状态轴划分过粗导致生态位宽度失去意义。我的处理是直接放弃该轴的宽度计算，而不是强行填一个0，否则数字看起来完整，实际毫无意义。

7.4 坑四：宽度和重叠被混为一谈

生态位宽度描述的是单个物种的资源利用范围，生态位重叠描述的是两个物种利用共享资源的程度。两个窄生态位物种如果有相同的偏好，它们之间的重叠可能很高；两个宽生态位物种如果偏好互补，重叠反而可能很低。论文里常见错误是拿宽度排序去推导种间竞争强度，逻辑上不成立。要分析重叠就另外算Planka指数或Morisita指数，写作上要把这两个概念清楚分开。

7.5 坑五：R包的输出结构看着像列表其实是矩阵

spaa包的niche.width返回值不算复杂，但初次使用很容易直接用res_levins["Ba"]取数，结果返回一个列表而不是向量，png画图时还会报错。我现在的固定操作是先跑str()和class()，确认结构后再取数。如果发现是矩阵，就写res_levins[, "Ba"]，如果是数据框，就用res_levins$Ba。不同spaa版本之间确实有这种细微差异，这也是我为什么在正式分析中保留手写版函数的原因之一。

7.6 坑六：抽样强度不均导致“假性广布”

某物种在海边、山地、城市绿地都有记录，看起来宽度极大，但仔细看数据，海边只调查了1次、城市绿地只有零星偶见，这样的“广布”是采样假象。我在一份大型监测数据里就见过这种情况，一个稀有物种因为零星记录覆盖了多个生境，Levins宽度竟排进前三。处理办法是：先看总记录量，低于某个阈值的物种直接不进宽度排序，或者用稀化法统一抽样强度。生态位宽度的比较前提是各物种的采样努力具有可比性，这一点无论如何强调都不过分。

我在实际项目中还有一个体会：生态位宽度很少单独支撑一篇论文，它更适合作为群落分析链条里的一环，和多样性指数、生态位重叠、排序分析放在一起，相互印证。计算本身半小时就能跑完，真正决定工作量的是数据质量控制和结果解释。你拿到的数据如果满足资源状态互斥、缺失值明确、抽样强度可比这三个条件，剩下的就是公式、函数和一张好图的事。