如何用Tidyverse替代嵌套循环实现更高效的Bootstrapping？

我编写了一段嵌套for循环代码，用于对长期部署的多台相机采集的物种数据执行Bootstrapping分析，需要针对n_cameras与days的每种组合运行100次。代码能正常输出预期结果，但运行效率偏低，希望借助Tidyverse工具实现相同逻辑，至少移除嵌套循环来提升运行效率。

现有代码如下：

library(tidyverse)
library(geosphere)
library(terra)

Sites <- c("Site1","Site2","Site3","Site4","Site5","Site6")
Dates <- seq(as.Date("2023-01-01"), as.Date("2023-01-15"), by="days")
Species <- c("Sp1","Sp2","Sp3")

lonlat <- data.frame(
  lon = c(15.0, 25.0, 12.0, 18.0, 19.0, 19.0),
  lat = c(22.0, 18.0, 25.0, 20.0, 18.0, 25.0),
  SiteID = Sites
)

df <- data.frame(
  SiteID = sample(Sites, 20, replace=TRUE),
  Date = sample(Dates, 20, replace=TRUE),
  Species = sample(Species, 20, replace=TRUE),
  Count = sample(0:3, 20, replace=TRUE)
) %>%
  merge(lonlat, by="SiteID")

# 原嵌套循环代码
count_df <- data.frame(matrix(nrow = length(Species)+2, ncol = 0))
count_df$Species <- c(Species,"days","n_cameras")

days <- c(10,15,20,30)
n_cameras <- c(2,3,4)
runs <- 10
day_counter <- 2

for(i in 1:runs){
  for(k in 1:length(n_cameras)) {
    for(j in 1:length(days)) {
      cam1 <- sample_n(lonlat, 1)
      lonlat[,4] <- unlist(transpose(as.list(distm(cam1[,1:2], lonlat[,1:2], fun = distHaversine))))
      closest <- lonlat[order(lonlat[,4])[1:n_cameras[k]],]
      subsample <- df[which(interaction(df$SiteID) %in% interaction(closest$SiteID)),]
      dates <- unique(subsample$Date)
      a <- sample(dates,1)
      b <- a + days[j]
      sample_data <- subsample[which(subsample$Date >= a & subsample$Date < b),]
      summary <- sample_data %>%
        group_by(Species) %>%
        summarise(count = length(Count > 1))
      combined <- count_df %>%
        left_join(summary, by = "Species") %>%
        mutate(count = coalesce(count, 0))
      combined$count[length(Species)+1] <- days[j]
      combined$count[length(Species)+2] <- n_cameras[k]
      count_df[,day_counter] <- combined$count
      day_counter <- day_counter + 1
    }
  }
}

优化思路 #

• 参数组合扁平化：用crossing生成所有测试参数的笛卡尔积（运行次数、相机数量、天数），彻底移除嵌套循环
• 封装单次逻辑：把单次Bootstrap的操作写成独立函数，逻辑清晰且便于批量调用
• 批量映射执行：用purrr::pmap_dfr对每个参数组合批量执行Bootstrap，自动整理结果为规范数据框
• 长格式存储优先：先以长格式保存结果（符合Tidyverse规范），按需转成原代码的宽格式

优化后代码 #

library(tidyverse)
library(geosphere)
library(terra)

# 生成模拟数据（与原代码一致）
Sites <- c("Site1","Site2","Site3","Site4","Site5","Site6")
Dates <- seq(as.Date("2023-01-01"), as.Date("2023-01-15"), by="days")
Species <- c("Sp1","Sp2","Sp3")

lonlat <- data.frame(
  lon = c(15.0, 25.0, 12.0, 18.0, 19.0, 19.0),
  lat = c(22.0, 18.0, 25.0, 20.0, 18.0, 25.0),
  SiteID = Sites
)

df <- data.frame(
  SiteID = sample(Sites, 20, replace=TRUE),
  Date = sample(Dates, 20, replace=TRUE),
  Species = sample(Species, 20, replace=TRUE),
  Count = sample(0:3, 20, replace=TRUE)
) %>%
  merge(lonlat, by="SiteID")

# 定义单次Bootstrap执行函数
bootstrap_single <- function(n_cam, day_length, df_data, lonlat_data, species_list) {
  # 随机选择起始相机
  cam1 <- sample_n(lonlat_data, 1)
  # 计算所有相机到起始相机的距离
  distances <- distm(cam1[,1:2], lonlat_data[,1:2], fun = distHaversine) %>% as.vector()
  # 筛选最近的n_cam个相机站点
  closest_sites <- lonlat_data %>% 
    mutate(dist = distances) %>% 
    arrange(dist) %>% 
    slice_head(n = n_cam) %>% 
    pull(SiteID)
  
  # 筛选对应站点的数据
  subsample <- df_data %>% filter(SiteID %in% closest_sites)
  
  # 处理无有效数据的边界情况
  if(nrow(subsample) == 0) {
    return(tibble(
      Species = c(species_list, "days", "n_cameras"),
      count = c(rep(0, length(species_list)), day_length, n_cam)
    ))
  }
  
  # 随机选择起始日期并筛选时间范围
  start_date <- sample(unique(subsample$Date), 1)
  end_date <- start_date + day_length
  sample_data <- subsample %>% filter(Date >= start_date, Date < end_date)
  
  # 统计物种记录数（与原代码逻辑一致：length(Count>1)等价于n()）
  species_counts <- sample_data %>%
    group_by(Species) %>%
    summarise(count = n(), .groups = "drop")
  
  # 补全所有未出现的物种计数为0
  full_counts <- tibble(Species = species_list) %>%
    left_join(species_counts, by = "Species") %>%
    mutate(count = coalesce(count, 0))
  
  # 添加days和n_cameras的标识记录
  bind_rows(
    full_counts,
    tibble(Species = "days", count = day_length),
    tibble(Species = "n_cameras", count = n_cam)
  )
}

# 设置分析参数
days <- c(10,15,20,30)
n_cameras <- c(2,3,4)
runs <- 10

# 生成所有参数组合
param_grid <- crossing(
  run_id = 1:runs,
  n_cam = n_cameras,
  day_length = days
)

# 批量执行Bootstrap并整理结果
bootstrap_results <- param_grid %>%
  mutate(
    result = pmap(list(n_cam, day_length), 
                  ~bootstrap_single(..1, ..2, df_data = df, lonlat_data = lonlat, species_list = Species))
  ) %>%
  unnest(result)

# 转成原代码的宽格式（可选）
wide_results <- bootstrap_results %>%
  pivot_wider(
    id_cols = Species,
    names_from = c(run_id, n_cam, day_length),
    values_from = count
  )

优化说明 #

• 效率提升：移除嵌套循环，用向量式操作和批量映射替代，减少全局变量修改的开销，Tidyverse内部优化进一步加速运行
• 可读性增强：逻辑拆分到独立函数中，参数清晰，便于后续维护和调试
• 鲁棒性提升：添加了无数据场景的处理，避免运行报错
• 灵活性更高：长格式结果适配更多后续分析需求，可一键转成原代码的宽格式

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