R中具有p值的相关矩阵答案

【问题标题】：correlation matrix with p-value in RR中具有p值的相关矩阵
【发布时间】：2018-05-20 13:31:10
【问题描述】：

假设我想要进行相关矩阵

library(dplyr)

data(iris)
iris %>% 
  select_if(is.numeric) %>%
  cor(y =iris$Petal.Width, method = "spearman") %>%  round(2)

现在我们看到了

              [,1]
Sepal.Length  0.83
Sepal.Width  -0.29
Petal.Length  0.94
Petal.Width   1.00

我希望统计显着相关性用 * 标记在哪里

*<0,05
**<0,01
*** <0,001

怎么办？

【问题讨论】：

你的意思是library(Hmisc);rcorr(as.matrix(iris[1:4]), iris$Petal.Width)$P
@D.Joe 请查看我的更新。我相信在我之前的帖子中我没有在正确的位置指定method = "spearman"，所以我得到了不同的P值。现在我已经解决了。

标签： r dataframe dplyr

【解决方案1】：

使用tidyverse 的解决方案。我们可以将数据框转换为长格式，使用nest创建列表列，然后使用map对每个子集执行cor.test。之后，map_dbl 可以通过指定名称"p.value" 来提取 P 值。 dat1 是最终输出。

library(tidyverse)

data(iris)
dat1 <- iris %>% 
  select_if(is.numeric) %>%
  gather(Column, Value, -Petal.Width) %>%
  group_by(Column) %>%
  nest() %>%
  mutate(Cor = map(data, ~cor.test(.x$Value, .x$Petal.Width, method = "spearman"))) %>%
  mutate(Estimate = round(map_dbl(Cor, "estimate"), 2), 
         P_Value = map_dbl(Cor, "p.value"))

dat1
# # A tibble: 3 x 5
#   Column       data               Cor         Estimate  P_Value
#   <chr>        <list>             <list>         <dbl>    <dbl>
# 1 Sepal.Length <tibble [150 x 2]> <S3: htest>    0.83  4.19e-40
# 2 Sepal.Width  <tibble [150 x 2]> <S3: htest>   -0.290 3.34e- 4
# 3 Petal.Length <tibble [150 x 2]> <S3: htest>    0.94  8.16e-70

如果您不需要列表列，可以使用select 将其删除。

dat1 %>% select(-data, -Cor)
# # A tibble: 3 x 3
#   Column       Estimate  P_Value
#   <chr>           <dbl>    <dbl>
# 1 Sepal.Length    0.83  4.19e-40
# 2 Sepal.Width    -0.290 3.34e- 4
# 3 Petal.Length    0.94  8.16e-70

现在我们可以使用mutate 和case_when 来添加标签以显示重要性。

dat2 <- dat1 %>%
  select(-data, -Cor) %>%
  mutate(Significance = case_when(
    P_Value < 0.001  ~ "*** <0,001",
    P_Value < 0.01   ~ "** <0,01",
    P_Value < 0.05   ~ "*<0,05",
    TRUE             ~ "Not Significant"
  ))
dat2
# # A tibble: 3 x 4
#   Column       Estimate  P_Value Significance
#   <chr>           <dbl>    <dbl> <chr>       
# 1 Sepal.Length    0.83  4.19e-40 *** <0,001  
# 2 Sepal.Width    -0.290 3.34e- 4 *** <0,001  
# 3 Petal.Length    0.94  8.16e-70 *** <0,001

【讨论】：

你也可以使用map2 map2(iris[1:3], iris["Petal.Width"], ~ {Cor <- cor.test(.x, .y, method = "spearman"); tibble(Estimate = Cor$estimate, PValue = Cor$p.value)}) %>% bind_rows(.id = 'Column')
@akrun 这是一个不错的解决方案。感谢分享。
当我得到 dat1 时，出现错误 dat1 %>% select(-data, -Cor) 选择错误 (., -data, -Cor) : 未使用的参数 (-data, -Cor ) 为什么它不起作用:)
@D.Joe 我不知道。你能确定你使用的是dplyr包中的select函数吗？

【解决方案2】：

这里有两个tidyverse 选项，它们都使用来自broom 的tidy。使用 tidy 将提取估计值和 p 值，因此您不必手动执行此操作。我为您想要显示的不同显着性水平制作了一个中断向量，因此您可以使用cut 轻松剪切和标记 p 值；将其保存在命名向量中也使其更具可重复性。

我第一次使用cor.test，它通过管道进入tidy.htest 方法。我第二次使用来自Hmisc 的rcorr，它通过管道进入tidy.rcorr 方法。

在第一种情况下，我将 gathered 数据框转换为长格式，以将每个度量与 Petal.Width 进行比较；在第二种情况下，需要一个矩阵，我使用了完整的数据集，然后过滤了包含Petal.Width 的任一列。

library(tidyverse)

sig_breaks <- c(zero = 0, "***" = 0.001, "**" = 0.01, "*" = 0.05, NS = Inf)

iris %>%
  as_tibble() %>%
  select_if(is.numeric) %>%
  gather(key = measure, value = value, -Petal.Width) %>%
  group_by(measure) %>%
  do(mtx = cor.test(.$value, .$Petal.Width, method = "spearman")) %>%
  broom::tidy(mtx) %>%
  mutate(stars = cut(p.value, breaks = sig_breaks, include.lowest = T, labels = names(sig_breaks)[2:5]))
#> # A tibble: 3 x 7
#> # Groups:   measure [3]
#>   measure      estimate statistic  p.value method        alternative stars
#>   <chr>           <dbl>     <dbl>    <dbl> <fct>         <fct>       <fct>
#> 1 Petal.Length    0.938    35061. 8.16e-70 Spearman's r… two.sided   ***  
#> 2 Sepal.Length    0.834    93208. 4.19e-40 Spearman's r… two.sided   ***  
#> 3 Sepal.Width    -0.289   725048. 3.34e- 4 Spearman's r… two.sided   ***

iris %>%
  select_if(is.numeric) %>%
  as.matrix() %>%
  Hmisc::rcorr(type = "spearman") %>%
  broom::tidy() %>%
  filter(column1 == "Petal.Width" | column2 == "Petal.Width") %>%
  mutate(stars = cut(p.value, breaks = sig_breaks, include.lowest = T, labels = names(sig_breaks)[2:5]))
#>        column1     column2   estimate   n      p.value stars
#> 1 Sepal.Length Petal.Width  0.8342888 150 0.0000000000   ***
#> 2  Sepal.Width Petal.Width -0.2890317 150 0.0003342981   ***
#> 3 Petal.Length Petal.Width  0.9376668 150 0.0000000000   ***

由reprex package (v0.2.0) 于 2018 年 5 月 20 日创建。

【讨论】：

【解决方案3】：

您可以根据自己的需要调整corstarsl()。

corFun <- function (x) {
  library(Hmisc)
  x <- as.matrix(x)
  R <- rcorr(x, type="spearman")$r
  p <- rcorr(x, type="spearman")$P
  stars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", 
                                             ifelse(p < 0.05, "* ", " ")))
  R <- format(round(cbind(rep(-1.11, ncol(x)), R), 2))[, -1]
  Rnew <- matrix(paste(R, stars, sep = ""), ncol = ncol(x))
  diag(Rnew) <- paste(diag(R), " ", sep = "")
  rownames(Rnew) <- colnames(x)
  colnames(Rnew) <- paste(colnames(x), "", sep = "")
  Rnew <- as.matrix(Rnew)
  Rnew <- as.data.frame(Rnew)
  return(Rnew)
}

收益

> data.frame(r=corFun(iris[, -5])[, 4])
                    r
Sepal.Length  0.83***
Sepal.Width  -0.29***
Petal.Length  0.94***
Petal.Width     1.00

【讨论】：