如何在R中对多组变量对的相关性进行置换检验？

我想用cor.test计算两个数据框中多组变量对的相关性，比较基于估计值或统计量得到的p值。真实数据情况：

• df1：20个样本、40个变量
• df2：20个样本、345个变量

现有双变量相关性置换示例多针对单变量对的for循环，我尝试用foreach对df2的每个y变量做1000次置换（测试阶段用10次+小数据集），通过置换统计量与原统计量对比计算p值，但代码中的置换逻辑未生效。测试代码如下：

df <- as.data.frame(matrix(round(runif(n=12*12, min=1, max=20), 0), ncol=12))
colnames(df) <- paste0("a", ncol(df))
df2 <- as.data.frame(matrix(round(runif(n=12*20, min=1, max=20), 0), ncol=20))
colnames(df2) <- paste("b", ncol(df2) )
comp <- expand.grid(colnames(df1), colnames(df2))


res <-  foreach(i=1:nrow(comp), .combine=rbind, .packages=c("magrittr", "dplyr")) %dopar% {
perm <- \(){
N <- 10
samp <- sample(nrow(df2))
cor.test(df[, comp[i,1]], df2[samp, comp[i,2]], method="spearman", exact=FALSE)[c("estimate")]
}
perm_stat <- replicate(N, perm())
orig_stat <- cor.test(df[, comp[i,1]], df2[, comp[i,2]], method="spearman", exact=FALSE)[c("estimate")] 
(sum(abs(perm_stat) >= abs(orig_stat)) + 1)/N 
   }

res2 <- as.data.frame(unlist(res)) %>% dplyr::mutate(Var1=comp[,1], Var2=comp[,2])

仅保留orig_stat时能正常计算相关系数，但置换逻辑失效。想问是否该采用foreach双重循环，或是mcmapply这类并行化Map方案？同时给出了双重循环示例代码：

##double loop format 
out <- foreach(j=1:ncol(df1), .combine = 'rbind', .packages=c("magrittr", "dplyr")) %:%
   foreach(i = 1:ncol(df2), .combine = 'c') %dopar% {
    a <- cor.test(df1[,j], df2[,i], method = "spearman")$p.value
  }

你的测试代码置换失效的核心原因是嵌套函数的变量作用域问题，同时存在测试数据生成的小错误（比如df2的列名全重复）。以下是修正后的可行方案：

修正后的单循环并行代码 #

先修复基础错误，再优化置换逻辑：

# 生成正确的测试数据
df1 <- as.data.frame(matrix(round(runif(n=12*12, min=1, max=20), 0), ncol=12))
colnames(df1) <- paste0("a", 1:ncol(df1))
df2 <- as.data.frame(matrix(round(runif(n=12*20, min=1, max=20), 0), ncol=20))
colnames(df2) <- paste0("b", 1:ncol(df2))
comp <- expand.grid(colnames(df1), colnames(df2))

# 加载并行工具包
library(foreach)
library(doParallel)
library(dplyr)

# 初始化并行集群
cl <- makeCluster(detectCores() - 1)
registerDoParallel(cl)

res <- foreach(i=1:nrow(comp), .combine=rbind, .packages=c("dplyr")) %dopar% {
  var1 <- df1[, comp[i,1]]
  var2 <- df2[, comp[i,2]]
  N <- 10
  
  # 计算原统计量
  orig_stat <- cor.test(var1, var2, method="spearman", exact=FALSE)$estimate
  
  # 直接在replicate内实现置换，避免函数作用域问题
  perm_stat <- replicate(N, {
    perm_var2 <- sample(var2)
    cor.test(var1, perm_var2, method="spearman", exact=FALSE)$estimate
  })
  
  # 计算置换p值（+1避免p=0，符合置换检验惯例）
  p_val <- (sum(abs(perm_stat) >= abs(orig_stat)) + 1) / (N + 1)
  data.frame(Var1=comp[i,1], Var2=comp[i,2], orig_r=orig_stat, perm_p=p_val)
}

# 关闭集群
stopCluster(cl)

# 查看结果
head(res)

关键修正点： #

• 移除嵌套的perm函数，直接在replicate内写置换逻辑，避免变量传递问题
• 修正df2的列名生成错误（原代码列名全为b 20，导致变量对重复）
• 直接对单个变量var2采样，比置换整个数据框行更高效
• 调整p值计算分母为N+1，符合置换检验的无偏估计标准

可选方案：双重循环并行 #

如果偏好更直观的双重循环逻辑，可使用foreach的嵌套循环（%:%）：

# 初始化集群
cl <- makeCluster(detectCores() - 1)
registerDoParallel(cl)

out <- foreach(j=1:ncol(df1), .combine='rbind', .packages=c("dplyr")) %:%
  foreach(i=1:ncol(df2), .combine='rbind') %dopar% {
    var1 <- df1[,j]
    var2 <- df2[,i]
    N <- 10
    
    orig_stat <- cor.test(var1, var2, method="spearman", exact=FALSE)$estimate
    perm_stat <- replicate(N, {
      cor.test(var1, sample(var2), method="spearman", exact=FALSE)$estimate
    })
    p_val <- (sum(abs(perm_stat) >= abs(orig_stat)) + 1) / (N + 1)
    
    data.frame(Var1=colnames(df1)[j], Var2=colnames(df2)[i], orig_r=orig_stat, perm_p=p_val)
  }

stopCluster(cl)

替代方案：mcmapply并行化映射 #

mcmapply是mapply的并行版本，适合处理多变量对的映射任务，代码更简洁：

library(parallel)

# 定义单个变量对的置换检验函数
perm_cor_test <- function(var1_name, var2_name) {
  var1 <- df1[, var1_name]
  var2 <- df2[, var2_name]
  N <- 10
  
  orig_stat <- cor.test(var1, var2, method="spearman", exact=FALSE)$estimate
  perm_stat <- replicate(N, {
    cor.test(var1, sample(var2), method="spearman", exact=FALSE)$estimate
  })
  p_val <- (sum(abs(perm_stat) >= abs(orig_stat)) + 1) / (N + 1)
  
  data.frame(Var1=var1_name, Var2=var2_name, orig_r=orig_stat, perm_p=p_val)
}

# 初始化并行集群
cl <- makeCluster(detectCores() - 1)
clusterExport(cl, c("df1", "df2"))  # 将数据传递到集群节点
res_mcmapply <- mcmapply(perm_cor_test, comp$Var1, comp$Var2, SIMPLIFY=FALSE, mc.cores=detectCores()-1)
res_mcmapply <- do.call(rbind, res_mcmapply)
stopCluster(cl)

方案选择建议： #

• 你的数据规模（40*345=13800对变量）下，三种方案均可行
• foreach单/双重循环更易调试，适合需要分步输出中间结果的场景
• mcmapply代码更简洁，适合逻辑清晰的映射任务

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