Obtaining a matrix of p-values of a pearson correlation matrix

1
votes

Dear StackOverflow Community:

I am trying to make a matrix of p-values that corresponds to a matrix I have obtained by obtaining correlation values

My data is the following (just doing 5 rows for simplicity, my real data is 3 columns for each data frame with 50 rows). 

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")

But I get this error:

Error in rcorr(t(FG_Smooth), t(FMG_Smooth), type = "pearson") : must have >4 observations

I only have 3 biological samples of each - so I am unable to use the rcorr command that has been suggested mulitple time in multiple posts. The rcorr command gives you 1) the matrix of correlations; and 2) p-values for the correlations.

So, to side-step this issue, I have ran the following: as suggested in other posts: 

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

This works and gives a matrix of all of my correlations.

My next step is to find the p-values associated with each value in the correlation matrix. The function cor.test only gives an overall p-value, which isn't what I need.

After perusing multiple posts - I ran across this one: rcorr() function for correlations

I followed directions to code given: 

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"

This only gives me 9 p-values and not a matrix of p-values that corresponds to each correlation value obtained with the cor command. For this example, it would be a 5x5 matrix of p-values, since the cor command results in a 5x5 matrix of correlation values.

Is there a diff. way to do this?

2
rcorrelationp-value

2 Answers

1
votes

Here's a tidyverse solution that creates all pairs of interest and then performs a cor.test for each pair and extracts the correlation value and the corresponding p value.

# 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 and v2 are the row names of your datasets that will create the pairs for the correlation tests, cor_test column has the correlation test object for each pair, cor_value has the extracted correlation coefficient and cor_p_value has the extracted p value.

If you save the above output as a data frame you can easily reshape. For example if you save it as d you can get a 5x5 data frame of p values like this:

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

An alternative version using broom package as well would be:

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

which gives you a tidy format of the correlation test object. You need to use columns estimate (correlation coefficient) and p.value.

0
votes

A couple of issues arise with your current setup:

  1. First your cor.test must used the transposed version with t():

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

  2. Second, cor.test returns a list of elements which you only need to extract the p.value item:

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

Therefore, consider binding the mapply results into a 5 X 5 matrix with needed colnames and rownames (i.e., dimnames):

P-Value matrix

# 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

Correlation matrix

And to be sure, if we use $estimate, the matrix build replicates cor() call:

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