如何生成含丰度的随机群落矩阵并批量计算多样性指数？

拥有不同样地的物种丰度数据，想要评估现有群落集合的多样性指数与随机生成群落的指数是否存在差异，目前面临两个核心问题：

1. 生成仅随机化物种组成但保留丰度数据（非仅存在/缺失型）的随机群落矩阵；
2. 批量计算这些随机群落的功能多样性指数（优先使用FD或mFD包），分类或系统发育多样性指数也可作为示例。

曾尝试使用vegan包的nullmodel和simulate函数，但生成的是存在/缺失型随机群落，对函数帮助文档理解困难，已尝试的代码如下：

library(vegan)
data("mite")
community <- mite
nm <- vegan::nullmodel(x = community, method = "swap") # 基于观测群落矩阵（35列×75行）生成空模型
sm <- simulate(nm, nsim = 999)

一、生成保留丰度的随机群落矩阵 #

你之前使用的swap方法默认针对存在/缺失矩阵，要保留丰度，需选择适配丰度数据的空模型方法：

1. 适配丰度的空模型方法 #

vegan的nullmodel支持多种适合丰度数据的随机化算法，常用选项：

• 'quasiswap'：针对丰度矩阵的交换算法，保留样地总丰度和物种出现频率
• 'swapcount'：基于丰度计数的交换，严格保留样地和物种的边际总和
• 'r2dtable'：生成服从行列边际总和约束的随机矩阵（适合严格控制总丰度的场景）

2. 修正后的代码示例 #

library(vegan)
data("mite")
community <- mite

# 选择适配丰度的空模型方法（此处以quasiswap为例）
nm_abund <- nullmodel(community, method = "quasiswap")
# 生成999个随机群落，返回结果为矩阵列表
sm_abund <- simulate(nm_abund, nsim = 999)

# 验证：查看第一个随机群落的丰度分布（非0/1格式）
head(sm_abund[[1]])

若需要将随机群落整理为三维数组（样地×物种×模拟次数），方便后续批量计算：

library(abind)
# 转换为三维数组
sm_array <- abind(sm_abund, along = 3)

二、批量计算多样性指数 #

以下分别展示功能多样性（FD包）和系统发育多样性（vegan包）的批量计算流程：

1. 功能多样性指数（FD包） #

首先需要准备物种功能性状数据，此处使用mite配套的性状数据作为示例：

library(FD)
data("mite.trait")

# 计算观测群落的功能多样性指数（此处以FDiv为例，可按需选择FRic、FEve等）
obs_fd <- dbFD(mite.trait, community, calc.FRic = TRUE, calc.FDiv = TRUE)
obs_fdiv <- obs_fd$FDiv

# 批量计算999个随机群落的FDiv指数
random_fdiv <- lapply(sm_abund, function(random_comm) {
  fd_result <- dbFD(mite.trait, random_comm, calc.FRic = TRUE, calc.FDiv = TRUE)
  return(fd_result$FDiv)
})

# 将结果转换为数据框，便于后续统计分析
random_fdiv_df <- do.call(rbind, random_fdiv)
colnames(random_fdiv_df) <- rownames(community)

2. 系统发育多样性指数（vegan包） #

若有物种的系统发育树，可使用vegan的pd函数计算系统发育多样性：

library(ape)
# 模拟系统发育树（实际使用时替换为你的真实物种树）
set.seed(123)
phy_tree <- rcoal(ncol(community))
phy_tree$tip.label <- colnames(community)

# 计算观测群落的系统发育多样性（PD）
obs_pd <- pd(community, phy_tree, include.root = FALSE)$PD

# 批量计算随机群落的PD指数
random_pd <- lapply(sm_abund, function(random_comm) {
  pd(random_comm, phy_tree, include.root = FALSE)$PD
})

random_pd_df <- do.call(rbind, random_pd)

3. 检验观测值与随机分布的差异 #

通过计算观测值在随机分布中的百分位数，判断是否显著偏离随机预期：

# 以FDiv为例，计算每个样地观测值的百分位数
fdiv_percentiles <- apply(random_fdiv_df, 2, function(random_vals) {
  ecdf(random_vals)(obs_fdiv)
})

# 输出结果：值接近0或1说明观测值显著偏离随机群落的指数分布
print(fdiv_percentiles)

内容的提问来源于stack exchange，提问作者Pedro Eça