如何在R中为Spearman相关矩阵构建BCa自助法置信区间？

使用boot包实现BCa自助法计算置信区间，结合rstatix获取匹配的p值，同时控制结果精度到四位小数，解决你遇到的三个问题：

1. 加载所需包 #

library(boot)
library(rstatix)
library(tidyverse)
library(ggcorrplot)

2. 准备示例数据 #

df_test <- data.frame(
    recall=c(45, 1, 32, 17, 79, 15, 75, 100, 43, 80, 74, 91, 60, 54, 67, 26, 97, 53, 51, 30),
    recog=c(90, 29, 93, 73, 34, 68, 78, 56, 92, 85, 35, 81, 7, 58, 4, 52, 82, 31, 6, 23),
    hits=c(77, 89, 44, 8, 70, 96, 76, 62, 95, 27, 49, 12, 16, 28, 83, 2, 36, 10, 61, 86),
    misses=c(59, 78, 14, 44, 86, 61, 80, 72, 25, 93, 5, 42, 64, 95, 73, 54, 7, 67, 11, 53),
    false_alarms=c(32, 34, 55, 41, 77, 76, 89, 36, 12, 100, 70, 62, 47, 81, 90, 63, 13, 83, 79, 38)
)

3. 定义自助统计量函数 #

该函数返回所有变量对的Spearman相关系数（转换为向量方便boot包处理）：

spearman_cor_fun <- function(data, indices) {
  boot_sample <- data[indices, ]
  cor_mat <- cor(boot_sample, method = "spearman", use = "pairwise.complete.obs")
  # 提取下三角的相关系数（排除对角线）
  cor_vec <- cor_mat[lower.tri(cor_mat)]
  return(cor_vec)
}

4. 运行BCa自助过程 #

指定10000次迭代：

boot_result <- boot(data = df_test, statistic = spearman_cor_fun, R = 10000)

# 提取所有变量对名称（下三角）
var_pairs <- expand.grid(colnames(df_test), colnames(df_test)) %>%
  filter(Var1 > Var2) %>%
  mutate(pair = paste(Var2, Var1, sep = "-")) %>%
  pull(pair)

5. 计算并整理BCa置信区间 #

对每个相关系数计算BCa CI，保留四位小数：

bca_ci_list <- list()

for (i in 1:length(var_pairs)) {
  ci <- boot.ci(boot_result, index = i, type = "bca")
  # 提取BCa上下限并保留四位小数
  bca_ci <- round(ci$bca[4:5], 4)
  bca_ci_list[[var_pairs[i]]] <- bca_ci
}

# 转换为对称矩阵形式
bca_ci_mat_lower <- matrix(NA, nrow = ncol(df_test), ncol = ncol(df_test))
rownames(bca_ci_mat_lower) <- colnames(df_test)
colnames(bca_ci_mat_lower) <- colnames(df_test)

for (pair in names(bca_ci_list)) {
  vars <- strsplit(pair, "-")[[1]]
  bca_ci_mat_lower[vars[2], vars[1]] <- paste0("[", bca_ci_list[[pair]][1], ", ", bca_ci_list[[pair]][2], "]")
}

# 生成完整的对称CI矩阵
bca_ci_mat <- t(bca_ci_mat_lower) + bca_ci_mat_lower
diag(bca_ci_mat) <- "[1.0000, 1.0000]"

6. 获取匹配的相关系数和p值（四位小数） #

# 计算Spearman相关矩阵
cor_mat <- round(cor(df_test, method = "spearman", use = "pairwise.complete.obs"), 4)

# 获取与rstatix匹配的p值矩阵
p_mat <- round(cor_pmat(df_test, method = "spearman"), 4)

7. 合并结果并展示 #

# 整合所有结果
result_list <- list(
  correlation_matrix = cor_mat,
  p_value_matrix = p_mat,
  bca_ci_matrix = bca_ci_mat
)

# 查看各结果
print(result_list$correlation_matrix)
print(result_list$p_value_matrix)
print(result_list$bca_ci_matrix)

8. 可视化（可选：结合CI与p值） #

# 生成带CI的标注文本
label_mat <- matrix(NA, nrow = ncol(df_test), ncol = ncol(df_test))
for (i in 1:ncol(df_test)) {
  for (j in 1:ncol(df_test)) {
    if (i > j) {
      label_mat[i, j] <- paste0(cor_mat[i,j], "\n", bca_ci_mat[i,j])
    } else if (i == j) {
      label_mat[i, j] <- "1.0000\n[1.0000, 1.0000]"
    }
  }
}

# 绘制带CI和p值的相关图
ggcorrplot(cor_mat,
           type = "lower", 
           method = "circle", 
           colors = c("turquoise3", "khaki1", "violetred3"), 
           p.mat = p_mat,
           lab = TRUE,
           lab_matrix = label_mat)

方案说明 #

1. p值匹配：直接使用rstatix::cor_pmat()生成p值，与你之前的结果完全一致，CI单独通过boot计算，避免不一致问题。
2. BCa自助法：通过boot.ci(..., type = "bca")明确指定BCa类型，满足需求。
3. 高精度小数：所有结果通过round(..., 4)强制保留四位小数，达到项目精度要求。

内容的提问来源于stack exchange，提问作者Linda Jasmine Hoffman