如何在R中按组进行相关性检验并将结果存储为列表？

我的数据如下：

x = data.frame(year = c(2019, 2019, 2020, 2020, 2021, 2021, 2022, 2022, 2023, 2023),
               group = rep(c("A", "B"), 5),
               y = runif(10))
x
#    year group          y
# 1  2019     A 0.26550866
# 2  2019     B 0.37212390
# 3  2020     A 0.57285336
# 4  2020     B 0.90820779
# 5  2021     A 0.20168193
# 6  2021     B 0.89838968
# 7  2022     A 0.94467527
# 8  2022     B 0.66079779
# 9  2023     A 0.62911404
# 10 2023     B 0.06178627

我需要对每个group组内的year和y变量做Pearson相关性检验，目前单独筛选分组能实现：

library(dplyr)
A = x %>% filter(group == "A")
B = x %>% filter(group == "B")

cor.test(A$year, A$y)
# 
#     Pearson's product-moment correlation
# 
# data:  A$year and A$y
# t = 0.67259, df = 3, p-value = 0.5494
# alternative hypothesis: true correlation is not equal to 0
# 95 percent confidence interval:
#  -0.7644073  0.9430671
# sample estimates:
#       cor 
# 0.3619872

cor.test(B$year, B$y)
# 
#     Pearson's product-moment correlation
# 
# data:  B$year and B$y
# t = 1.9909, df = 3, p-value = 0.1406
# alternative hypothesis: true correlation is not equal to 0
# 95 percent confidence interval:
#  -0.3822520  0.9826436
# sample estimates:
#       cor 
# 0.7544519 

但尝试批量处理多分组时两次都失败了：

失败尝试1：用dplyr的group_map #

x %>% group_by(group) %>% 
  group_map(~cor.test(x$year, x$y))

# [[1]]
# 
#     Pearson's product-moment correlation
# 
# data:  x$year and x$y
# t = 0.15447, df = 8, p-value = 0.8811
# alternative hypothesis: true correlation is not equal to 0
# 95 percent confidence interval:
#  -0.5955411  0.6614485
# sample estimates:
#        cor 
# 0.05453354 
# 
# 
# [[2]]
# 
#     Pearson's product-moment correlation
# 
# data:  x$year and x$y
# t = 0.15447, df = 8, p-value = 0.8811
# alternative hypothesis: true correlation is not equal to 0
# 95 percent confidence interval:
#  -0.5955411  0.6614485
# sample estimates:
#        cor 
# 0.05453354

结果不符合预期，因为代码错误引用了原数据集x，而非分组后的子数据框。

失败尝试2：用plyr的daply #

library(plyr)
daply(x, .(group), function(y) cor.test(y$year, y$y))
# 执行报错

解决方案 #

方法1：修正dplyr的group_map用法 #

在group_map的匿名函数中，用.指代当前分组的子数据框，而非原数据集x：

x %>% 
  group_by(group) %>% 
  group_map(~cor.test(.$year, .$y))

# [[1]]
# 
#     Pearson's product-moment correlation
# 
# data:  .$year and .$y
# t = 0.67259, df = 3, p-value = 0.5494
# alternative hypothesis: true correlation is not equal to 0
# 95 percent confidence interval:
#  -0.7644073  0.9430671
# sample estimates:
#       cor 
# 0.3619872
# 
# 
# [[2]]
# 
#     Pearson's product-moment correlation
# 
# data:  .$year and .$y
# t = 1.9909, df = 3, p-value = 0.1406
# alternative hypothesis: true correlation is not equal to 0
# 95 percent confidence interval:
#  -0.3822520  0.9826436
# sample estimates:
#       cor 
# 0.7544519

方法2：用dplyr的summarise存储结果为列表列 #

如果需要保留分组信息，可将结果存入列表列，方便后续调用：

result_df <- x %>% 
  group_by(group) %>% 
  summarise(cor_result = list(cor.test(year, y)))

# 查看分组A的结果
result_df$cor_result[[1]]

方法3：用purrr的split+map #

先按group拆分数据集，再逐个处理：

library(purrr)
x_split <- split(x, x$group)
map(x_split, ~cor.test(.$year, .$y))

以上三种方法都能得到每个分组的相关性检验结果列表。

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