How to run multiple semi-partial correlations in R?-CodePudding

I am a beginner in R. What I want to do is compute semi partial correlations in R, while controlling for one variable.

For example, in the case of the mtcars data set, I would like to compute all the correlations between cyl, disp, hp, drat, wt, qsec, vs, am, gear, carb, while controlling for mpg.

From what I observed, most packages, including ppcor, do something like this:

pcor.test(a, b, c)

Which means computing the correlation between a and b while controlling for c. But I do not know of a way to extend that to multiple variables while controlling for one. I am not sure how to extend this syntax to achieve that. In pseudocode, what I would want is something like:

compute_correlations(cyl, disp, hp, drat, wt, qsec, vs, am, gear, carb; control for = mpg)

And then visualize the result in a heatmap, with labeling that says how I computed all these correlations while controlling for 'mpg'. Is there any way to do this efficiently in R? I guess what I could do is make all these correlations myself for each pair of variables, but that would take a long time and it would be inefficient.

This is quite easily achievable in SPSS but I do not know how to do this in R.

CodePudding user response：

The process to do this is indeed not evident, and perhaps this is an opportunity to create a package that makes it easier. In the meantime, here is one solution where we iterate through each combination of variables:

library(ppcor)
library(purrr)

mtcars_tidy <- tibble::rownames_to_column(mtcars)

vars_to_correlate <- names(mtcars_tidy)[3:12] |> 
  combn(2) |> 
  t() |> 
  as.data.frame()

map2_df(vars_to_correlate[[1]], vars_to_correlate[[2]], 
        ~ pcor.test(mtcars_tidy[[.x]], mtcars_tidy[[.y]], mtcars_tidy[["mpg"]])) |> 
  cbind(vars_to_correlate)

      estimate        p.value  statistic  n gp  Method   V1   V2
1   0.64737094 0.000082740213  4.5740136 32  1 pearson  cyl disp
2   0.51835526 0.002814957618  3.2641973 32  1 pearson  cyl   hp
3  -0.31184550 0.087670622303 -1.7674786 32  1 pearson  cyl drat
4   0.16570850 0.372986653785  0.9048777 32  1 pearson  cyl   wt
5  -0.49337691 0.004796854818 -3.0545737 32  1 pearson  cyl qsec
6  -0.62604621 0.000165211632 -4.3234453 32  1 pearson  cyl   vs
7  -0.02735369 0.883868197138 -0.1473593 32  1 pearson  cyl   am
8  -0.18172321 0.327880737991 -0.9951794 32  1 pearson  cyl gear
9   0.13169303 0.480070640823  0.7154196 32  1 pearson  cyl carb
10  0.39777260 0.026688732321  2.3347229 32  1 pearson disp   hp
11 -0.34200622 0.059673996837 -1.9599493 32  1 pearson disp drat
12  0.57833847 0.000654820766  3.8176756 32  1 pearson disp   wt
13 -0.16358622 0.379224930550 -0.8929679 32  1 pearson disp qsec
14 -0.37198547 0.039339008725 -2.1580693 32  1 pearson disp   vs
15 -0.19507995 0.292950712364 -1.0711167 32  1 pearson disp   am
16 -0.31902184 0.080242339460 -1.8127037 32  1 pearson disp gear
17 -0.16247217 0.382523780650 -0.8867212 32  1 pearson disp carb
18  0.17318251 0.351501338785  0.9469246 32  1 pearson   hp drat
19 -0.04690019 0.802172856767 -0.2528435 32  1 pearson   hp   wt
20 -0.66932190 0.000038342372 -4.8513350 32  1 pearson   hp qsec
21 -0.44054279 0.013124294295 -2.6426556 32  1 pearson   hp   vs
22  0.44076669 0.013072569401  2.6443225 32  1 pearson   hp   am
23  0.44676093 0.011749935870  2.6891769 32  1 pearson   hp gear
24  0.61230870 0.000251252635  4.1706361 32  1 pearson   hp carb
25 -0.33357553 0.066672764314 -1.9054998 32  1 pearson drat   wt
26 -0.29177519 0.111237981706 -1.6427381 32  1 pearson drat qsec
27 -0.02200590 0.906461717887 -0.1185341 32  1 pearson drat   vs
28  0.51916508 0.002764864921  3.2711766 32  1 pearson drat   am
29  0.58000604 0.000626263781  3.8342517 32  1 pearson drat gear
30  0.46560757 0.008299545471  2.8332183 32  1 pearson drat carb
31  0.41764131 0.019400101391  2.4752780 32  1 pearson   wt qsec
32  0.05714957 0.760085155921  0.3082637 32  1 pearson   wt   vs
33 -0.43250238 0.015098604771 -2.5831987 32  1 pearson   wt   am
34 -0.38196965 0.033966232663 -2.2257362 32  1 pearson   wt gear
35 -0.12149573 0.514995012825 -0.6591576 32  1 pearson   wt carb
36  0.68705027 0.000019643524  5.0919744 32  1 pearson qsec   vs
37 -0.66197095 0.000049951046 -4.7560806 32  1 pearson qsec   am
38 -0.51946280 0.002746643801 -3.2737454 32  1 pearson qsec gear
39 -0.56154329 0.001012819623 -3.6546170 32  1 pearson qsec carb
40 -0.38439791 0.032753726959 -2.2423299 32  1 pearson   vs   am
41 -0.17215886 0.354399233697 -0.9411560 32  1 pearson   vs gear
42 -0.32655948 0.072971566382 -1.8605799 32  1 pearson   vs carb
43  0.72095777 0.000004762368  5.6025889 32  1 pearson   am gear
44  0.58105260 0.000608906812  3.8446919 32  1 pearson   am carb
45  0.73589464 0.000002387013  5.8528090 32  1 pearson gear carb

EDIT: To create a heatmap:

cors_df$V1 <- factor(cors_df$V1, levels = names(mtcars_tidy)[3:12])
cors_df$V2 <- factor(cors_df$V2, levels = names(mtcars_tidy)[3:12])

library(ggplot2)

ggplot(cors_df, aes(x = V1, y = V2, fill = estimate, label = round(p.value, 4)))   geom_tile()   geom_text()