获得 pearson 相关矩阵的 p 值矩阵

Question

尊敬的 * 社区：

我正在尝试制作一个 p 值矩阵，该矩阵对应于我通过获得相关值获得的矩阵

我的数据如下（为简单起见只做 5 行，我的真实数据是 3 列，每个数据框有 50 行）。

FG_Smooth <- data.frame(FS_1 = c(0.43, 0.33, 3.47, 5.26, 1.09), FS2 = c(0.01, 0.02, 6.86, 3.27, 0.86), FS_3 = c(0.07, 0.36, 1.91, 5.61, 0.84), row.names = c("Group_3", "Thermo", "Embryophyta", "Flavo", "Cyclo")) 
FMG_Smooth <- data.frame(GS_1 = c(1.13, 1.20, 0.52, 2.81, 0.70), GS_2 = c(1.18, 1.7, 0.42, 2.93, 0.78), GS_3 = c(1.17, 1.11, 0.60, 3.10, 0.87), row.names = c("Proline", "Trigonelline", "L-Lysine", "Nioctine", "Caffeate"))                         

library(Hmisc)
rcorr(t(FG_Smooth), t(FMG_Smooth), type = "pearson")

但我得到这个错误：

rcorr(t(FG_Smooth), t(FMG_Smooth), type = "pearson") 中的错误： 必须有 >4 个观察结果

我每个人只有 3 个生物样本 - 所以我无法使用在多个帖子中多次建议的 rcorr 命令。 rcorr 命令为您提供 1) 相关矩阵； 2) 相关性的 p 值。

所以，为了回避这个问题，我运行了以下内容：正如其他帖子中所建议的那样：

library(stats)
cor(t(FG_Smooth), t(FMG_Smooth), method = "pearson")

这可行，并给出了我所有相关性的矩阵。

我的下一步是找到与相关矩阵中的每个值相关联的 p 值。函数cor.test 只给出了一个整体 p 值，这不是我需要的。

在阅读了多篇帖子后 - 我遇到了这个： rcorr() function for correlations

我按照给定代码的说明进行操作：

tblcols <- expand.grid(1:ncol(FG_Smooth), 1:ncol(FMG_Smooth))
cfunc <- function(var1, var2) {
  cor.test(FG_Smooth[,var1], FMG_Smooth[,var2], method="pearson")
}

res <- mapply(function(a,b) {
  cfunc(var1 = a, var2 = b)
}, tblcols$Var1, tblcols$Var2)

head(res)
[,1]        [,2]        [,3]        [,4]        [,5]        [,6]        [,7]       
statistic   1.324125    -0.1022017  2.422883    0.9131595   -0.3509424  1.734178    1.53494    
parameter   3           3           3           3           3           3           3          
p.value     0.2773076   0.9250449   0.09392613  0.4284906   0.74883     0.1812997   0.2223626  
estimate    0.6073388   -0.05890371 0.8135079   0.4663678   -0.1985814  0.7075406   0.663238   
null.value  0           0           0           0           0           0           0          
alternative "two.sided" "two.sided" "two.sided" "two.sided" "two.sided" "two.sided" "two.sided"
            [,8]         [,9]       
statistic   -0.009291327 2.880821   
parameter   3            3          
p.value     0.99317      0.06348644 
estimate    -0.005364273 0.8570256  
null.value  0            0          
alternative "two.sided"  "two.sided"

这只会给我 9 个 p 值，而不是对应于使用 cor 命令获得的每个相关值的 p 值矩阵。对于此示例，它将是一个 5x5 的 p 值矩阵，因为 cor 命令会生成一个 5x5 的相关值矩阵。

有区别吗？怎么做？

Answer 1

这是一个 tidyverse 解决方案，它创建所有感兴趣的对，然后为每对执行 cor.test 并提取相关值和相应的 p 值。

# example data
FG_Smooth <- data.frame(FS_1 = c(0.43, 0.33, 3.47, 5.26, 1.09), FS2 = c(0.01, 0.02, 6.86, 3.27, 0.86), FS_3 = c(0.07, 0.36, 1.91, 5.61, 0.84), row.names = c("Group_3", "Thermo", "Embryophyta", "Flavo", "Cyclo")) 
FMG_Smooth <- data.frame(GS_1 = c(1.13, 1.20, 0.52, 2.81, 0.70), GS_2 = c(1.18, 1.7, 0.42, 2.93, 0.78), GS_3 = c(1.17, 1.11, 0.60, 3.10, 0.87), row.names = c("Proline", "Trigonelline", "L-Lysine", "Nioctine", "Caffeate"))                         

library(tidyverse)

expand.grid(v1 = row.names(FG_Smooth),                                # create combinations of names
            v2 = row.names(FMG_Smooth)) %>%
  tbl_df() %>%                                                        # use for visualisation purpose
  mutate(cor_test = map2(v1, v2, ~cor.test(unlist(FG_Smooth[.x,]),    # perform the correlation test for each pair and store it
                                           unlist(FMG_Smooth[.y,]))), 
         cor_value = map_dbl(cor_test, "estimate"),                   # get the correlation value from the test
         cor_p_value = map_dbl(cor_test, "p.value"))                  # get the p value from the test

# # A tibble: 25 x 5
#   v1          v2           cor_test    cor_value cor_p_value
#   <fct>       <fct>        <list>          <dbl>       <dbl>
# 1 Group_3     Proline      <S3: htest>    -0.998     0.0367 
# 2 Thermo      Proline      <S3: htest>    -0.592     0.596  
# 3 Embryophyta Proline      <S3: htest>     0.390     0.745  
# 4 Flavo       Proline      <S3: htest>    -0.544     0.634  
# 5 Cyclo       Proline      <S3: htest>    -0.966     0.167  
# 6 Group_3     Trigonelline <S3: htest>    -0.492     0.673  
# 7 Thermo      Trigonelline <S3: htest>    -0.998     0.0396 
# 8 Embryophyta Trigonelline <S3: htest>     0.985     0.109  
# 9 Flavo       Trigonelline <S3: htest>    -1.000     0.00188
#10 Cyclo       Trigonelline <S3: htest>    -0.305     0.803  
# # ... with 15 more rows

v1 和 v2 是数据集的行名，它们将为相关测试创建对，cor_test 列包含每对的相关测试对象，cor_value 具有提取的相关系数，@ 987654328@有提取出来的p值。

如果您将上述输出保存为数据框，您可以轻松地进行整形。例如，如果您将其保存为d，您可以获得一个 5x5 的 p 值数据框，如下所示：

d %>%
  select(v1, v2, cor_p_value) %>%
  spread(v2, cor_p_value)

# # A tibble: 5 x 6
#   v1          Caffeate `L-Lysine` Nioctine Proline Trigonelline
#   <fct>          <dbl>      <dbl>    <dbl>   <dbl>        <dbl>
# 1 Cyclo          0.309     0.995     0.351  0.167       0.803  
# 2 Embryophyta    0.779     0.0931    0.737  0.745       0.109  
# 3 Flavo          0.890     0.204     0.848  0.634       0.00188
# 4 Group_3        0.439     0.875     0.481  0.0367      0.673  
# 5 Thermo         0.928     0.242     0.886  0.596       0.0396 

使用broom 包的替代版本是：

library(tidyverse)
library(broom)

expand.grid(v1 = row.names(FG_Smooth),                     
            v2 = row.names(FMG_Smooth)) %>%
  tbl_df() %>%
  mutate(cor_test = map2(v1, v2, ~tidy(cor.test(unlist(FG_Smooth[.x,]),    
                                                unlist(FMG_Smooth[.y,]))))) %>%
  unnest()

# # A tibble: 25 x 8
#   v1          v2           estimate statistic p.value parameter method                               alternative
#   <fct>       <fct>           <dbl>     <dbl>   <dbl>     <int> <chr>                                <chr>      
# 1 Group_3     Proline        -0.998   -17.3   0.0367          1 Pearson's product-moment correlation two.sided  
# 2 Thermo      Proline        -0.592    -0.735 0.596           1 Pearson's product-moment correlation two.sided  
# 3 Embryophyta Proline         0.390     0.423 0.745           1 Pearson's product-moment correlation two.sided  
# 4 Flavo       Proline        -0.544    -0.648 0.634           1 Pearson's product-moment correlation two.sided  
# 5 Cyclo       Proline        -0.966    -3.73  0.167           1 Pearson's product-moment correlation two.sided  
# 6 Group_3     Trigonelline   -0.492    -0.565 0.673           1 Pearson's product-moment correlation two.sided  
# 7 Thermo      Trigonelline   -0.998   -16.0   0.0396          1 Pearson's product-moment correlation two.sided  
# 8 Embryophyta Trigonelline    0.985     5.78  0.109           1 Pearson's product-moment correlation two.sided  
# 9 Flavo       Trigonelline   -1.000  -339.    0.00188         1 Pearson's product-moment correlation two.sided  
#10 Cyclo       Trigonelline   -0.305    -0.320 0.803           1 Pearson's product-moment correlation two.sided  
# # ... with 15 more rows

它为您提供了相关测试对象的tidy 格式。您需要使用列estimate（相关系数）和p.value。

Answer 2

您当前的设置存在几个问题：

1. 首先您的cor.test 必须使用带有t() 的转置版本：

   cor.test(t(FG_Smooth)[,var1], t(FMG_Smooth)[,var2], method="pearson")
   

2. 其次，cor.test返回一个元素列表，你只需要提取p.value项即可：

   cor.test(t(FG_Smooth)[,var1], t(FMG_Smooth)[,var2], method="pearson")$p.value
   

因此，请考虑将 mapply 结果绑定到一个 5 X 5 矩阵中，其中包含所需的 colnames 和 rownames（即 dimnames）：

P 值矩阵

# COMBINATION PAIRS
tblcols <- expand.grid(1:ncol(t(FG_Smooth)), 1:ncol(t(FMG_Smooth)))

# COR TEST FUNCTION
cfunc <- function(var1, var2) {
  cor.test(t(FG_Smooth)[,var1], t(FMG_Smooth)[,var2], method="pearson")$p.value
}

# P-VALUE MATRIX BUILD
matrix(mapply(cfunc, tblcols$Var1, tblcols$Var2),
       ncol = ncol(t(FG_Smooth)), nrow = ncol(t(FMG_Smooth)),
       dimnames = list(colnames(t(FG_Smooth)), colnames(t(FMG_Smooth))))

#               Proline Trigonelline   L-Lysine  Nioctine  Caffeate
# Group_3     0.0367145  0.672775387 0.87489349 0.4808196 0.4392690
# Thermo      0.5964129  0.039648033 0.24176614 0.8860530 0.9276037
# Embryophyta 0.7450881  0.109027230 0.09309087 0.7373778 0.7789284
# Flavo       0.6341827  0.001878145 0.20399625 0.8482831 0.8898338
# Cyclo       0.1669023  0.802963162 0.99491874 0.3506318 0.3090812

相关矩阵

可以肯定的是，如果我们使用 $estimate，矩阵构建会复制 cor() 调用：

t_FG <- t(FG_Smooth)
t_FMG <- t(FMG_Smooth)

cfunc <- function(var1, var2) {
  cor.test(t_FG[,var1], t_FMG[,var2], method="pearson")$estimate
}

# COR MATRIX BUILD
m <- matrix(mapply(cfunc, tblcols$Var1, tblcols$Var2),
            ncol = ncol(t_FG), nrow = ncol(t_FMG),
            dimnames = list(colnames(t_FG), colnames(t_FMG)))

cor(t(FG_Smooth), t(FMG_Smooth), method = "pearson")
#                Proline Trigonelline     L-Lysine   Nioctine   Caffeate
# Group_3     -0.9983375   -0.4916671  0.195254411 -0.7280867 -0.7712447
# Thermo      -0.5923344   -0.9980613  0.928751644  0.1780333  0.1134750
# Embryophyta  0.3898002    0.9853709 -0.989327898 -0.4009247 -0.3403212
# Flavo       -0.5435196   -0.9999956  0.949098001  0.2360668  0.1721863
# Cyclo       -0.9658300   -0.3045869 -0.007981547 -0.8521212 -0.8844400

m
#                Proline Trigonelline     L-Lysine   Nioctine   Caffeate
# Group_3     -0.9983375   -0.4916671  0.195254411 -0.7280867 -0.7712447
# Thermo      -0.5923344   -0.9980613  0.928751644  0.1780333  0.1134750
# Embryophyta  0.3898002    0.9853709 -0.989327898 -0.4009247 -0.3403212
# Flavo       -0.5435196   -0.9999956  0.949098001  0.2360668  0.1721863
# Cyclo       -0.9658300   -0.3045869 -0.007981547 -0.8521212 -0.8844400