如何 运行 R 中的多个半偏相关?

How to run multiple semi-partial correlations in R?

我是 R 的初学者。我想做的是在 R 中计算半偏相关,同时控制一个变量。

例如,对于 mtcars 数据集,我想计算 cyl、disp、hp、drat、wt、qsec、vs、am、gear、carb 之间的所有相关性,同时控制 mpg .

根据我的观察,大多数软件包,包括 ppcor,都会做这样的事情:

pcor.test(a, b, c)

这意味着在控制c的同时计算a和b之间的相关性。 但我不知道有什么方法可以在控制一个变量的同时将其扩展到多个变量。我不确定如何扩展此语法来实现它。在伪代码中,我想要的是:

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

然后在热图中可视化结果,标签说明我如何在控制 'mpg' 的同时计算所有这些相关性。有没有办法在 R 中有效地做到这一点?我想我能做的是自己为每对变量建立所有这些相关性,但这会花费很长时间,而且效率很低。

这在 SPSS 中很容易实现,但我不知道如何在 R 中实现。

这样做的过程确实不明显,也许这是创建一个使它更容易的包的机会。同时,这是我们迭代每个变量组合的解决方案:

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

编辑:创建热图:

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