如何自动找到多个变量的偏相关答案

【问题标题】：How to find the partial correlation for severals variables automatically如何自动找到多个变量的偏相关
【发布时间】：2021-06-30 07:11:54
【问题描述】：

我有几个变量。我想在给定第三个变量的情况下找到2 变量之间的部分相关性。但是，我需要计算所有可能的排列。也就是说，假设我有X1、X2、X2、X4。然后，我需要计算偏相关如下：

   X1,X2;X3
   X1,X2;X4,
   X1,X3;X2,
   X1,X2;X4,
   X1,X3;X4
   X1,X4;X2

等等。

我用来查找偏相关的函数是R中的pcor.test()表单包ppcor。

我的数据是：

structure(c(16.09, 15.74, 15.73, 16.1, 16.58, 16.54, 32.84, 32.06, 
32.1, 32.77, 32.36, 32.1, 33.64, 32.68, 33.7, 34.91, 36.57, 37.31, 
20.27, 19.79, 20.33, 20.5, 21.69, 22.72, 786, 765, 753, 757.5, 
748.5, 735, 45.49, 44.27, 44.42, 44.92, 46.36, 47.35, 29.46, 
28.5, 29.18, 29.66, 30.94, 30.06, 49.6, 48.76, 49.05, 49.08, 
48.95, 48.7, 49.1, 47.8, 48.35, 48.85, 48.92, 49.1, 46.28, 45.31, 
45.21, 48.04, 48.6, 49.15, 141, 136.5, 136.2, 137.5, 137.9, 136.5, 
1489, 1448, 1442, 1427, 1443, 1429, 354.2, 354, 362, 376.9, 377.6, 
381.3, 38.6, 37.49, 37.94, 38.85, 40.21, 41.11, 364, 358.1, 357.1, 
356.4, 359.3, 358.5, 253.5, 250, 250, 254, 258.4, 253.8, 145.3, 
143.9, 144, 145.5, 145.7, 146.8, 48.32, 47.03, 49.02, 50, 50.9, 
51.4, 10.55, 10.3, 10.91, 10.93, 9.75, 10.26, 166.7, 159.5, 163.7, 
165.6, 173, 174.3), .Dim = c(6L, 20L), .Dimnames = list(c("2011-09-09", 
"2011-09-12", "2011-09-13", "2011-09-14", "2011-09-15", "2011-09-16"
), c("X1", "X2", "X3", "X4", "X5", "X6", "X7", 
"X8", "X9", "X10", "X11", "X12", "X13", "X14", "X15", 
"X16", "X17", "X18", "X19", "X20")))

为了找到偏相关，我使用了这种方法：

pcor.test(dat[,1],dat[,2],dat[,1], method="kendall")

【问题讨论】：

标签： r combinations

【解决方案1】：

强制as.data.frame后直接在combn中应用ppcor::pcor.test。提供一系列不错的结果。

library(ppcor)
res <- combn(as.data.frame(dat), 3, function(x) pcor.test(x[1], x[2], x[3]))
res[,,1]  ## first result
# [[1]]
# [1] 0.5488149
# 
# [[2]]
# [1] 0.338091
# 
# [[3]]
# [1] 1.137128
# 
# [[4]]
# [1] 6
# 
# [[5]]
# [1] 1
# 
# [[6]]
# [1] "pearson"

要知道之后你可以做什么：

res1 <- unlist(combn(colnames(dat), 3, function(x) 
  setNames(
    list(do.call(pcor.test, unname(as.data.frame(dat[, x])))), 
    paste(x, collapse='.')),
  simplify=F), recursive=F)

res1$X1.X2.X3
#    estimate  p.value statistic n gp  Method
# 1 0.5488149 0.338091  1.137128 6  1 pearson

或者这个：

res2 <- do.call(rbind, combn(colnames(dat), 3, function(x) 
  cbind(do.call(pcor.test, unname(as.data.frame(dat[, x]))), 
        comb=paste0(x, c('+', ':', ''), collapse='')), simplify=F))

head(res2)
#     estimate   p.value  statistic n gp  Method     comb
# 1  0.5488149 0.3380910  1.1371279 6  1 pearson X1+X2:X3
# 2  0.7378542 0.1546276  1.8934466 6  1 pearson X1+X2:X4
# 3  0.6127243 0.2718743  1.3428747 6  1 pearson X1+X2:X5
# 4  0.4509408 0.4459447  0.8750759 6  1 pearson X1+X2:X6
# 5 -0.1104844 0.8596136 -0.1925434 6  1 pearson X1+X2:X7
# 6  0.4441672 0.4536571  0.8586700 6  1 pearson X1+X2:X8

编辑

您也可以选择method='kendall'。

res3 <- do.call(rbind, combn(colnames(dat), 3, function(x) 
  cbind(
    tryCatch(do.call(pcor.test, c(unname(as.data.frame(dat[, x])), 
                                  method='kendall')),
             error=function(e) {
               setNames(as.data.frame(t(rep(NA, 6))), 
                        c("estimate", "p.value", "statistic", "n", "gp", "Method"))
               }), 
        comb=paste0(x, c('+', ':', ''), collapse='')), simplify=F))

res3[304:308,]
#       estimate   p.value  statistic  n gp  Method       comb
# 304  0.1247236 0.7599785  0.3055091  6  1 kendall X2+X14:X15
# 305 -0.1604088 0.6943787 -0.3929198  6  1 kendall X2+X14:X16
# 306         NA        NA         NA NA NA    <NA> X2+X14:X17
# 307  0.2766417 0.4980057  0.6776309  6  1 kendall X2+X14:X18
# 308  0.1706640 0.6759180  0.4180398  6  1 kendall X2+X14:X19

【讨论】：

@AnoushiravanR 问题可能是可能出现的奇点。在这种情况下，您可以使用try.Catch 方法，请参阅编辑。
@AnoushiravanR 这与您的原始问题无关，我建议您在尝试捕获时提出一个新问题:)
@AnoushiravanR 你知道，我不使用 tidyverse。
我明白，我必须解决这个问题，直到我做对了。不过谢谢你的好建议。
@AnoushiravanR 哦对不起xD

【解决方案2】：

将我的朋友awesome answer 简化一点，在tidyverse 中给出的结果与@jay.sf's answer 的结果相同

library(tidyverse)
library(ppcor)

combn(colnames(df), 3) %>% 
  as.data.frame() %>% 
  as.list() %>%
  map_dfr(~ {dd <- df[,colnames(df) %in% .x]; pcor.test(dd[,1], dd[,2], dd[,3], method = 'pearson')}  %>%
            mutate(cols = paste0(paste(colnames(dd)[1:2], collapse = '+'),';', colnames(dd)[3] ))) %>%
  head()
#>     estimate   p.value  statistic n gp  Method     cols
#> 1  0.5488149 0.3380910  1.1371279 6  1 pearson X1+X2;X3
#> 2  0.7378542 0.1546276  1.8934466 6  1 pearson X1+X2;X4
#> 3  0.6127243 0.2718743  1.3428747 6  1 pearson X1+X2;X5
#> 4  0.4509408 0.4459447  0.8750759 6  1 pearson X1+X2;X6
#> 5 -0.1104844 0.8596136 -0.1925434 6  1 pearson X1+X2;X7
#> 6  0.4441672 0.4536571  0.8586700 6  1 pearson X1+X2;X8

对于kendall，可以这样尝试

combn(colnames(df), 3) %>% 
  as.data.frame() %>% 
  as.list() %>%
  map_dfr(~ {dd <- df[,colnames(df) %in% .x]; tryCatch(error = function(cond){setNames(as.data.frame(t(rep(NA, 6))),
                                                                                       c("estimate", "p.value", "statistic", "n", "gp", "Method"))}, 
                                                       pcor.test(dd[,1], dd[,2], dd[,3], method = 'kendall'))}  %>%
            mutate(cols = paste0(paste(colnames(dd)[1:2], collapse = '+'),';', colnames(dd)[3] ))) %>%
  head()

1  0.17251639 0.6726038  0.42257713 6  1 kendall X1+X2;X3
2  0.17251639 0.6726038  0.42257713 6  1 kendall X1+X2;X4
3  0.25383654 0.5340931  0.62177000 6  1 kendall X1+X2;X5
4 -0.03589791 0.9299311 -0.08793156 6  1 kendall X1+X2;X6
5 -0.21107063 0.6051455 -0.51701534 6  1 kendall X1+X2;X7
6  0.38924947 0.3403557  0.95346259 6  1 kendall X1+X2;X8

检查

combn(colnames(df), 3) %>% 
  as.data.frame() %>% 
  as.list() %>%
  map_dfr(~ {dd <- df[,colnames(df) %in% .x]; tryCatch(error = function(cond){setNames(as.data.frame(t(rep(NA, 6))),
                                                                                       c("estimate", "p.value", "statistic", "n", "gp", "Method"))}, 
                                                       pcor.test(dd[,1], dd[,2], dd[,3], method = 'kendall'))}  %>%
            mutate(cols = paste0(paste(colnames(dd)[1:2], collapse = '+'),';', colnames(dd)[3] ))) %>%
  {.[304:308,]}

      estimate   p.value  statistic  n gp  Method       cols
304  0.1247236 0.7599785  0.3055091  6  1 kendall X2+X14;X15
305 -0.1604088 0.6943787 -0.3929198  6  1 kendall X2+X14;X16
306         NA        NA         NA NA NA    <NA> X2+X14;X17
307  0.2766417 0.4980057  0.6776309  6  1 kendall X2+X14;X18
308  0.1706640 0.6759180  0.4180398  6  1 kendall X2+X14;X19

【讨论】：

【解决方案3】：

更新的解决方案 感谢亲爱的@Jay.sf 提出的关于在发生错误时使用tryCatch 控制输出的宝贵建议，以及我亲爱的朋友@AnilGoyal 的大力帮助，他真是太棒了。

您也可以使用以下解决方案，但是当方法设置为kendall 时，会引发以下错误：

solve.default(cvx) 中的错误：系统在计算上是奇异的：倒数条件数 = 2.77556e-17

我做了一些研究，发现这通常发生在矩阵是奇异且不可逆的时候。

解决方案 正如我之前提到的，我们可以将函数包装在tryCatch 中，以便在出现错误警告或消息时自定义输出。在这里，我们明确要求函数产生NA 值，而不是在发生错误时停止函数的执行：

library(purrr)
library(magrittr)

combn(colnames(df), 3) %>% 
  as.data.frame() %>% 
  as.list() %>%
  map(~ df[names(df) %in% .x]) %>%
  set_names({.} %>% map(~ names(.x) %>% paste(collapse = " "))) %>%
  map(~ tryCatch(pcor.test(.x[1], .x[2], .x[3], method = "kendall"),
             error = function(cond) {
               setNames(as.data.frame(t(rep(NA, 6))),
                        c("estimate", "p.value", "statistic", "n", "gp", "Method"))
             })) %>%
  imap_dfr(~ .x %>% 
             mutate(Vars = .y)) %>%
  slice_head(n = 5)

     estimate   p.value   statistic n gp  Method     Vars
1  0.17251639 0.6726038  0.42257713 6  1 kendall X1 X2 X3
2  0.17251639 0.6726038  0.42257713 6  1 kendall X1 X2 X4
3  0.25383654 0.5340931  0.62177000 6  1 kendall X1 X2 X5
4 -0.03589791 0.9299311 -0.08793156 6  1 kendall X1 X2 X6
5 -0.21107063 0.6051455 -0.51701534 6  1 kendall X1 X2 X7

【讨论】：

【解决方案4】：

主要问题是生成所有可能的组合，这可以通过combn 完成

tmp=combn(colnames(df),3)

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]  [,9]  [,10] ...
[1,] "X1" "X1" "X1" "X1" "X1" "X1" "X1" "X1"  "X1"  "X1"  ...
[2,] "X2" "X2" "X2" "X2" "X2" "X2" "X2" "X2"  "X2"  "X2"  ...
[3,] "X3" "X4" "X5" "X6" "X7" "X8" "X9" "X10" "X11" "X12" ...

您可以从这里循环（或应用）该矩阵并执行计算。

【讨论】：