如何在R中实现两数据框相关系数的Bootstrap并行计算

核心思路 #

• 先计算所有列对的原始Spearman相关统计量
• 定义Bootstrap抽样函数：对两个数据框同步进行有放回行抽样，计算该抽样下所有列对的相关系数
• 并行执行指定次数的Bootstrap迭代，收集所有迭代的统计结果
• 基于Bootstrap结果计算p值：统计每个列对中，Bootstrap统计量绝对值超过原始统计量的次数占总迭代次数的比例

基础准备代码 #

# 加载依赖包
library(dplyr)

# 生成示例数据（固定随机种子保证可复现）
set.seed(123)
D1 <- data.frame(matrix(runif(10*10, 0, 2), ncol=10)) 
D2 <- data.frame(matrix(runif(10*16, 0, 2), ncol=16)) 

colnames(D1) <- paste0("a", 1:ncol(D1))
colnames(D2) <- paste0("b", 1:ncol(D2))

# 生成所有列对组合
compare <- expand.grid(colnames(D1), colnames(D2))

# 计算原始相关统计量（仅保留后续需要的估计值和统计量）
original_cor <- Map(function(x,y) {
  ct <- cor.test(D1[, x], D2[, y], method="spearman")
  c(estimate = ct$estimate, statistic = ct$statistic)
}, compare$Var1, compare$Var2) %>% 
  sapply(unlist) %>% 
  t() %>% 
  as.data.frame() %>% 
  mutate(Var1 = compare$Var1, Var2 = compare$Var2)

# Bootstrap迭代次数
R <- 100

Unix HPC平台（用parallel::mclapply实现并行） #

library(parallel)

# 定义单次Bootstrap迭代函数
bootstrap_iter <- function(seed) {
  set.seed(seed) # 为每个迭代设置独立种子，保证结果可复现
  # 同步抽样：两个数据框使用相同行索引，保留样本配对关系
  sample_idx <- sample(nrow(D1), replace = TRUE)
  D1_boot <- D1[sample_idx, ]
  D2_boot <- D2[sample_idx, ]
  
  # 计算该抽样下所有列对的统计量
  boot_cor <- Map(function(x,y) {
    ct <- cor.test(D1_boot[, x], D2_boot[, y], method="spearman")
    ct$statistic
  }, compare$Var1, compare$Var2) %>% 
    unlist()
  
  return(boot_cor)
}

# 设置并行核数（根据HPC资源调整，建议留1核给系统）
n_cores <- detectCores() - 1

# 并行执行Bootstrap
boot_results <- mclapply(1:R, bootstrap_iter, mc.cores = n_cores)

# 转换结果为矩阵（每行对应一次迭代的所有列对统计量）
boot_stats_matrix <- do.call(rbind, boot_results)

# 计算Bootstrap p值
original_stats <- original_cor$statistic
bootstrap_pvals <- apply(boot_stats_matrix, 2, function(boot_stats) {
  mean(abs(boot_stats) >= abs(original_stats[col(boot_stats_matrix)[1]]))
})

# 合并p值到原始结果
final_result <- original_cor %>% 
  mutate(bootstrap_p = bootstrap_pvals)

head(final_result)

Windows平台（用parallel::clusterMap实现并行） #

library(parallel)

# 创建并行集群
n_cores <- detectCores() - 1
cl <- makeCluster(n_cores)

# 导出集群节点需要的对象
clusterExport(cl, c("D1", "D2", "compare"))

# 定义单次Bootstrap迭代函数
bootstrap_iter_win <- function(seed) {
  set.seed(seed)
  sample_idx <- sample(nrow(D1), replace = TRUE)
  D1_boot <- D1[sample_idx, ]
  D2_boot <- D2[sample_idx, ]
  
  boot_cor <- Map(function(x,y) {
    ct <- cor.test(D1_boot[, x], D2_boot[, y], method="spearman")
    ct$statistic
  }, compare$Var1, compare$Var2) %>% 
    unlist()
  
  return(boot_cor)
}

# 并行执行Bootstrap
boot_results_win <- clusterMap(cl, bootstrap_iter_win, seed = 1:R)

# 关闭集群释放资源
stopCluster(cl)

# 转换结果并计算p值
boot_stats_matrix_win <- do.call(rbind, boot_results_win)
bootstrap_pvals_win <- apply(boot_stats_matrix_win, 2, function(boot_stats) {
  mean(abs(boot_stats) >= abs(original_stats[col(boot_stats_matrix_win)[1]]))
})

# 合并结果
final_result_win <- original_cor %>% 
  mutate(bootstrap_p = bootstrap_pvals_win)

head(final_result_win)

补充说明 #

• 同步抽样：D1和D2的行对应同一观测样本，必须使用相同的抽样索引，避免破坏配对关系
• 可复现性：为每个迭代设置独立随机种子，解决并行环境下随机数冲突问题
• 资源适配：Unix平台mclapply是fork式并行，更轻量；Windows平台需要显式创建集群，用clusterMap执行任务
• 性能优化：超大型数据集可拆分列对分组并行，减少单任务计算量

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