使用 group_by 创建相关性和 p 值的数据框，然后在 R 中绘制误差条

Question

我想绘制数据集中几个因素之间的相关性。如果可能，我想尝试为这些绘制的值添加误差条或晶须。在计算值之前，我首先想根据其中一列中的值对它们进行分组。如果可能的话，我想使用 tidyverse 解决方案。我可以使用 cor() 实现其中一半，但我不知道如何添加包含 p 值的附加列。

我认为 iris 数据集很好地展示了我想做的事情。实际数据使用沿 x 轴的时间序列。我指定 spearman 是因为这是我分析中使用的相关性，而不是因为它是 iris 数据集中的正确选择。我看到其他一些 posts 建议使用 cor.test 并从中提取值，但我不确定如何将其应用回条形图以用作误差线。下面是创建基本条形图的代码。

Edit 我已经将示例从使用 mtcars 数据集更改为 iris 数据集因为我认为它更好地反映了我的数据。虽然 jay.sf 对问题的最初回答适用于 mtcars 集并且非常感谢，但它不适用于我的数据集并且 iris 集抛出了与我相同的错误有。此外，我没有在原文中说明这一点，但 tidyverse 解决方案更可取，但不是必需的。

我认为我正在寻找的答案可能包含在这里，但我仍在努力找出细节：https://dominicroye.github.io/en/2019/tidy-correlation-tests-in-r/。

iristest <- iris %>%
  group_by(Species) %>%
  summarise(COR = cor(Sepal.Length,Sepal.Width, method = "spearman", use="complete.obs"))

ggplot(data = iristest) +
  aes(x = Species, y = COR) +
  geom_bar(stat = "identity") +
  theme_minimal()

原样，iristest 提供此输出：

    Species     COR
1   setosa      0.7553375
2   versicolor  0.5176060
3   virginica   0.4265165

我认为理想情况下，我希望输出在 COR 列之后添加 p 值。

    Species     COR          p-value
1   setosa      0.7553375    ###
2   versicolor  0.5176060    ###
3   virginica   0.4265165    ###

Answer 1

cor.test 生成一个列表，其中实际上存储了您需要的所有内容。因此，只需编写一个获取所需值的函数即可。我们可以在这里使用 by ，它会产生一个列表，我们可以 rbind 得到一个具有完美行名称的矩阵以进行绘图。列表数据帧的 rbind 需要 do.call。

res <- do.call(rbind, by(iris, iris$Species, function(x) { 
  rr <- with(x, cor.test(Sepal.Length, Sepal.Width, method="pearson")) 
  return(c(rr$estimate, CI=rr$conf.int)) 
})) 
#                  cor       CI1       CI2
# setosa     0.7425467 0.5851391 0.8460314
# versicolor 0.5259107 0.2900175 0.7015599
# virginica  0.4572278 0.2049657 0.6525292

请注意，method="spearman" 不适用于 iris 等相关数据，因此我在这里使用了 "pearson"。

为了绘制数据，我推荐 R 附带的 barplot。我们存储条形位置 b <- 并将它们用作 arrows 的 x 坐标。对于 y 坐标，我们从矩阵中获取值。

b <- barplot(res[,1], ylim=c(0, range(res)[2]*1.1), 
             main="My Plot", xlab="cyl", ylab="Cor. Sepal.Length ~ Sepal.Width")
arrows(b, res[,2], b, res[,3], code=3, angle=90, length=.1)
abline(h=0)
box()

Answer 2

主要使用 tidyverse...

这是与 Spearman 完成的关联：

library(tidyverse)
library(RVAideMemoire)

iristest <-  iris %>%
+ group_by(Species) %>%
+ group_modify(~ glance(spearman.ci(.x$Sepal.Width, .x$Sepal.Length))


iristest
# A tibble: 3 x 5
# Groups:   Species [3]
  Species    estimate conf.low.Inf conf.high.Sup method                     
  <fct>         <dbl>        <dbl>         <dbl> <chr>                      
1 setosa        0.755        0.599         0.857 Spearman's rank correlation
2 versicolor    0.518        0.251         0.724 Spearman's rank correlation
3 virginica     0.427        0.131         0.653 Spearman's rank correlation

使用 ggplot...

ggplot(iristest, aes(x = Species, y = estimate)) 
+ geom_bar(stat="identity") 
+ geom_errorbar(aes(ymin=conf.low.Inf, ymax=conf.high.Sup), width=.2, position=position_dodge(.9))

Answer 3

这是一个实现所要求的版本。分解成步骤，它比上面的例子稍微长一些。此版本仅使用基础 R，但某些人可能会感兴趣。

# Just extract the columns used in your question
data = iris[, c("Sepal.Length", "Sepal.Width", "Species")]

# Group the data by species
grouped.data = by(data, (data$Species), list)
# Run the function 'cor.test' (from stats) over the data from each species
cor.results = lapply(grouped.data, function(x) cor.test(x$Sepal.Length, x$Sepal.Width, method = "spearman", exact = FALSE) )
# Extract the rho and p-value
rho = sapply(cor.results, "[[", "estimate"))
p = sapply(cor.results, "[[", "p.value")
# Bundle the results into a data.frame (or whatever data structure you prefer)
data.frame(Species = names(cor.results), COR = rho, `p-value` = p, row.names = NULL)

     Species       COR      p.value
1     setosa 0.7553375 2.316710e-10
2 versicolor 0.5176060 1.183863e-04
3  virginica 0.4265165 2.010675e-03

[请参阅 ?cor.test 中有关使用这些数据所必需的 exact = FALSE 的注释。]

使用 group_by 创建相关性和 p 值的数据框，然后在 R 中绘制误差条

Create dataframe of correlation and p values using group_by then plot with error bars in R

r

tidyverse

correlation

errorbar

standard-error