How to find the partial correlation for severals variables automatically

Admin 1
Viewed 207

I have several variables. I would like to find the partial correlation between 2 variables given the third one. However, I need to calculate all possible permutations. That is, suppose I have X1, X2, X2, X4. Then, I need to compute the partial correlation as following:

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

and so on.

The function I used to find the partial correlation is pcor.test() form package ppcor in R.

My data is:

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

To find the partial correlation, I used this method:

pcor.test(dat[,1],dat[,2],dat[,1], method="kendall")
combinations r
4 Answers

After coercing as.data.frame apply ppcor::pcor.test directly in combn. Giving a nice array of results.

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"

To know which is what afterwards you could do this:

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

Or this:

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

Edit

You can also choose 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
user_65740380

Updated Solution Thanks to a valuable suggestion by dear @Jay.sf on using tryCatch to control the output while an error occurs and also great help by my dear friend @AnilGoyal who is so brilliant to be true.

You can also use the following solution, however when the method is set to kendall it throws the following error:

Error in solve.default(cvx) : system is computationally singular: reciprocal condition number = 2.77556e-17

I did some research and realized this usually happens when the matrix is singular and not invertible.

Solution As I mentioned earlier we can wrap our function in tryCatch in order to customize the output while an error a warning or a message occurs. Here we explicitly ask the function to produce NA values instead of bringing the execution of the function to a halt while an error occurs:

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
user_143145200

Simplifying my friends awesome answer a bit, in tidyverse which is giving same results as those by @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

For kendall try it like this

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

Checking it

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
user_28848590

The main issue is generating all the possible combinations, which can be done with 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" ...

from here you can loop (or apply) over this matrix and perform your calculation.

user_29749510
Related