Получение матрицы p-значений корреляционной матрицы Пирсона - PullRequest
Новая
       75

Получение матрицы p-значений корреляционной матрицы Пирсона

0 голосов
/ 19 ноября 2018

Уважаемое сообщество StackOverflow:

Я пытаюсь создать матрицу значений 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) р-значения для корреляций.

Итак, чтобы обойти эту проблему, я выполнил следующее: как предлагалось в других сообщениях: 

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

Это работает и дает матрицу всех моих корреляций.

Мой следующий шаг - найти p-значения, связанные с каждым значением в корреляционной матрице. Функция cor.test дает только общее значение p, а это не то, что мне нужно.

После просмотра нескольких постов - я наткнулся на это: Функция rcorr () для корреляций

Я следовал инструкциям к указанному коду: 

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-значений, а не матрицу p-значений, которая соответствует каждому значению корреляции, полученному с помощью команды cor. В этом примере это будет матрица p-значений 5x5, поскольку команда cor дает матрицу значений корреляции 5x5.

Есть ли разница? способ сделать это?

Ответы [ 2 ]

0 голосов
/ 19 ноября 2018

При текущей настройке возникает пара проблем:

  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() 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
0 голосов
/ 19 ноября 2018

Вот решение 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 имеет извлеченный коэффициент корреляции, а cor_p_value имеет извлеченное значение 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.

...