使用dplyr、tidyverse和broom生成带P值的相关矩阵

你提到的问题很常见：cor()函数只能返回相关系数矩阵，但不会给出显著性检验的P值。要同时拿到相关系数和P值，我们可以结合dplyr、purrr（tidyverse的一部分）和broom包来批量处理所有变量对，得到结构化的整洁结果。

方法一：用cor.test() + broom::tidy() 批量处理 #

这种方法灵活可控，能直接得到每对变量的相关系数和P值：

首先加载所需包，修正你的示例数据（统一样本量，避免不必要的缺失）：

set.seed(1164)
library(tidyverse)
library(broom)

# 修正示例数据（统一10个样本，标准差改为0.5避免误解）
ds <- data.frame(
  id = 1:10,
  a = rnorm(10, 2, 1),
  b = rnorm(10, 3, 2),
  c = rnorm(10, 1, 0.5)
)

接下来生成所有非重复的变量对，批量运行Spearman检验并提取结果：

# 提取要分析的变量名
target_vars <- ds %>% select(a, b, c) %>% colnames()

# 生成所有非重复的两两变量组合（排除自身相关和反向重复）
var_pairs <- expand.grid(var1 = target_vars, var2 = target_vars, stringsAsFactors = FALSE) %>%
  filter(var1 < var2)

# 批量运行检验并整理成整洁数据框
cor_results <- var_pairs %>%
  # 对每对变量运行Spearman相关检验
  mutate(test = map2(var1, var2, ~ cor.test(ds[[.x]], ds[[.y]], method = "spearman", use = "pairwise.complete.obs"))) %>%
  # 用broom把检验结果转换成数据框格式
  mutate(tidy_results = map(test, tidy)) %>%
  # 展开嵌套的数据框
  unnest(tidy_results) %>%
  # 保留核心列
  select(var1, var2, estimate, p.value, method)

# 查看结果
cor_results

运行后你会得到结构化的输出，包含变量对、相关系数(estimate)、P值(p.value)：

# A tibble: 3 × 4
  var1  var2  estimate p.value method                    
  <chr> <chr>    <dbl>   <dbl> <chr>                     
1 a     b       -0.245   0.473 Spearman's rank correlation rho
2 a     c        0.455   0.181 Spearman's rank correlation rho
3 b     c       -0.127   0.727 Spearman's rank correlation rho

方法二：用psych::corr.test()快速得到矩阵结果 #

如果你更习惯矩阵形式的输出，可以用psych包的corr.test()函数，它直接返回相关系数和P值矩阵，再用dplyr整理：

library(psych)

# 运行Spearman相关检验，同时得到相关系数和P值
cor_test_output <- corr.test(ds %>% select(a, b, c), method = "spearman", use = "pairwise.complete.obs")

# 整理相关系数矩阵为整洁数据框
cor_matrix_tidy <- cor_test_output$r %>%
  as.data.frame() %>%
  rownames_to_column("var1") %>%
  pivot_longer(cols = -var1, names_to = "var2", values_to = "correlation")

# 整理P值矩阵为整洁数据框
p_matrix_tidy <- cor_test_output$p %>%
  as.data.frame() %>%
  rownames_to_column("var1") %>%
  pivot_longer(cols = -var1, names_to = "var2", values_to = "p_value")

# 合并两个结果
combined_results <- cor_matrix_tidy %>%
  inner_join(p_matrix_tidy, by = c("var1", "var2")) %>%
  filter(var1 != var2)  # 排除变量自身的相关

combined_results

适配你的实际数据集 #

把代码替换成你的agreg_base_tipo_a数据集即可，还能快速筛选特定变量对的结果：

# 针对你的实际数据
target_vars <- agreg_base_tipo_a %>% select(S2.RT, BIS_total, IDATE, BAI, ASRS_total) %>% colnames()

var_pairs <- expand.grid(var1 = target_vars, var2 = target_vars, stringsAsFactors = FALSE) %>%
  filter(var1 < var2)

cor_results <- var_pairs %>%
  mutate(test = map2(var1, var2, ~ cor.test(agreg_base_tipo_a[[.x]], agreg_base_tipo_a[[.y]], method = "spearman", use = "pairwise.complete.obs"))) %>%
  mutate(tidy_results = map(test, tidy)) %>%
  unnest(tidy_results) %>%
  select(var1, var2, estimate, p.value)

# 快速查看特定变量对的结果，比如S2.RT和BIS_total
cor_results %>% filter(var1 == "S2.RT" & var2 == "BIS_total")

为什么你之前的代码得不到P值？ #

因为cor()函数仅计算相关系数，不执行显著性检验。要获取P值，必须使用cor.test()（针对单个变量对）或corr.test()（批量）这类能输出完整检验结果的函数，再用broom把非结构化的检验结果转换成适合dplyr处理的整洁数据框。

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