基于geiger包sim.bd函数的1000次模拟：按时间求物种丰富度均值

核心思路 #

每次sim.bd模拟返回包含time和n的时间序列数据，我们需要先收集1000次模拟的完整记录，再按time分组计算对应n的平均值。关键是确保每次模拟返回全时间点数据，再通过分组聚合得到均值。

完整代码示例 #

1. 加载包与设置参数

library(geiger)
library(dplyr) # 可选，用于更简洁的分组操作

# 自定义模拟参数
lambda <- 0.2  # 物种形成率
mu <- 0.08     # 灭绝率
max_time <- 15 # 总模拟时长
sim_times <- 1000 # 模拟次数

2. 批量模拟并存储结果

# 初始化空列表存储每轮模拟结果
sim_list <- vector("list", length = sim_times)

for (i in 1:sim_times) {
  # 必须设置return.all=TRUE，才能获取每个时间点的物种数变化
  single_sim <- sim.bd(n = 1, lambda = lambda, mu = mu, time = max_time, return.all = TRUE)
  # 转换为数据框存入列表
  sim_list[[i]] <- as.data.frame(single_sim)
}

3. 合并数据并计算时间点均值

方式一：使用dplyr（推荐，代码更直观）

# 合并所有模拟数据
all_sim_data <- bind_rows(sim_list)

# 按time分组，计算每组n的平均值
mean_n_time <- all_sim_data %>%
  group_by(time) %>%
  summarise(mean_species = mean(n, na.rm = TRUE))

# 查看结果
print(mean_n_time)

方式二：使用基础R

# 合并所有模拟数据
all_sim_base <- do.call(rbind, sim_list)

# 按time分组计算均值
mean_n_base <- aggregate(n ~ time, data = all_sim_base, FUN = mean, na.rm = TRUE)

# 查看结果
print(mean_n_base)

关键注意事项 #

• 必须添加return.all=TRUE参数：默认sim.bd只返回最终时间点的物种数，开启这个参数才能得到完整的时间序列。
• 用列表存储模拟结果：避免在循环中直接合并数据框，一次性合并列表的效率远高于循环内反复rbind。
• 替换错误的length()：你之前用length()得到的是该时间点的模拟次数总和，正确计算平均值要用mean()函数。

内容的提问来源于stack exchange，提问作者Georgia H