如何引导两个数据帧之间的相关性？

Question

我无法找到两个数据帧之间自举相关性的示例。我看过的关键帖子是 (1) https://stats.stackexchange.com/questions/20701/computing-p-value-using-bootstrap-with-r/ (2) 如何使用应用于大矩阵的向量化函数来引导相关性？

这里的NB/ p值需要通过计算绝对bootstrap检验统计量超出理论值多少倍的比例来获得，

幸运的是，我最近还遇到了一个使用 Map 来应用两个数据框之间的相关性的示例。 如何使用lapply代替嵌套的for循环来获取两个数据帧之间的相关性？

我有大型数据集，将在基于 UNIX 的 HPC 中运行，还有一个 Windows 操作系统选项，用于在 R 中对较小数据集运行计算。

   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))

  need_modify <- Map(function(x,y) cor.test(D1[, x], D2[, y], method="spearman"), compare$Var1,    compare$Var2) %>% 
  lapply(`[`, c('estimate', 'statistic', 'p.value')) %>% 
  sapply(unlist) %>%  t() %>% as.data.frame() %>% mutate(Var1=compare$Var1, Var2=compare$Var2)

  boot_df <- function(x) x[sample(nrow(x), rep = T), ]

 #number of bootstraps 
 R <- 100

我一直坚持修改上述内容，以便它能够使用基于 Unix 的操作系统（mclapply 或 mcMap）的并行化以及 Windows 的单独并行化（clusterMap 或 future_mapply）成功运行。

感谢任何正确方向的指点或其他地方的示例。

Answer 1

从技术上讲，您可以通过替换重新采样 B 次，提取

"statistic"

，

stopifnot(identical(dim(D1)[1], dim(D2)[1]))  ## check identical nrows
r.names <- Reduce(\(...) paste(..., sep=':'), Cmp)  ## store Cmp names

B <- 999L
set.seed(42)  ## for sake of reproducibility
res_boot <- parallel::mcMap(function(x, y) {
  boot <- function() {  ## boot fun
    smp <- sample.int(dim(D1)[1], dim(D1)[1],  ## resample w/ replacement
                      replace=TRUE)
    unlist(cor.test(D1[smp, x], D2[smp, y], method="spearman", exact=FALSE)[
      c('statistic')])  ## extract statistic
  }
  replicate(B, boot())  ## replicate `boot()` fun `B` times
}, Cmp$Var1, Cmp$Var2)  |> setNames(r.names)

lapply(res_boot, head) |> head(3)
# $`a1:b1`
# statistic.S statistic.S statistic.S statistic.S statistic.S statistic.S 
#    244.4444    277.0370    111.3750    242.4074    214.6012    194.3558 
# 
# $`a2:b1`
# statistic.S statistic.S statistic.S statistic.S statistic.S statistic.S 
#    97.54717   149.81595    37.59494   106.51899   199.62963   279.23077 
# 
# $`a3:b1`
# statistic.S statistic.S statistic.S statistic.S statistic.S statistic.S 
#    271.1465    145.9615    223.4161    261.1656    195.5556    187.6471 

并将超出实际统计量的次数加起来除以 B（通过将

+1

添加到分子和分母来进行偏差校正）。

res <- parallel::mcMap(\(x, y) cor.test(D1[, x], D2[, y], meth='s'), Cmp$Var1, Cmp$Var2) |>
  lapply(`[`, c('estimate', 'statistic', 'p.value')) |> 
  sapply(unlist) |> t() |> `rownames<-`(r.names)
head(res, 3)
#      estimate.rho statistic.S   p.value
# a1:b1  -0.27272727         210 0.4482722
# a2:b1   0.05454545         156 0.8916388
# a3:b1  -0.40606061         232 0.2473708

p_values <- mapply(\(x, y) (sum(unlist(x) > y) + 1)/(B + 1), 
       res_boot, res[, 'statistic.S']) |>
  setNames(r.names)

head(p_values)
# a1:b1 a2:b1 a3:b1 a4:b1 a5:b1 a6:b1 
# 0.516 0.446 0.506 0.510 0.493 0.519

请注意，答案重点在于如何技术上做到这一点。如果您确实希望它是合理的，请阅读戴维森和欣克利 (1997) 或咨询统计学家。

---

数据：

set.seed(42)
D1 <- data.frame(matrix(runif(10*10, 0, 2), ncol=10)) |> 
  setNames(paste0("a", seq_len(10)))
D2 <- data.frame(matrix(runif(10*16, 0, 2), ncol=16)) |> 
  setNames(paste0("a", seq_len(16)))
Cmp <- expand.grid(colnames(D1), colnames(D2))