Using a test for Multiple Corrections (e.g., BH) in a Partial Correlations with Multiple being Partialed Out

[nviet]

Hey all!

I'll keep this concise. I want to estimate a partial correlation using multiple variables (See Pcorr be- low) but include a correction for multiple comparisons like seen in the corr.test below ("adjust="BH""). Is there a package out there that allows this? I tried {pcorr} (obviously) and {corrr} which I think is just for visuals - but found nada. The code below will generate everything but I also left you all with a reprex too. Thank you all so much!

library(tidyverse)
library(ppcor)

#gendata

set.seed(123)
num_rows <- 20
num_cols <- 10
data <- matrix(rnorm(num_rows * num_cols), nrow = num_rows)
df <- as.data.frame(data) |>
set_names(c("Col1", "Col2", "Col3", "Col4", "Col5", "Col6", "Col7", "Col8", "Col9", "Col10"))

#Correlation

corr.test(df, y = NULL, use = "pairwise",method="pearson",adjust="BH",
alpha=.05,ci=TRUE,minlength=5,normal=TRUE)

#Pcorr

pcor.test(df$Col1, df$Col2, list(df$Col3, df$Col4, df$Col5, df$Col6, df$Col7, df$Col8, df$Col9, df$Col10), method = "pearson")

#Can I throw a multiple comparison correction in here ^ like "adjust="BH" in the #Correlation above?

tibble::tribble(
             ~Col1........Col2........Col3........Col4........Col5.......Col6........Col7........Col8........Col9.......Col10,
  "1   2.19881035 -0.57397348 -0.78862197 -0.52111732 -1.66747510 -0.7152422  0.23743027  0.62418747 -0.20529926  0.03455107",
  "2   1.31241298  0.61798582 -0.50219872 -0.48987045  0.73649596 -0.7526890  1.21810861  0.95900538  0.65119328  0.19023032",
  "3  -0.26514506  1.10984814  1.49606067  0.04715443  0.38602657 -0.9385387 -1.33877429  1.67105483  0.27376649  0.17472640",
  "4   0.54319406  0.70758835 -1.13730362  1.30019868 -0.26565163 -1.0525133  0.66082030  0.05601673  1.02467323 -1.05501704",
  "5  -0.41433995 -0.36365730 -0.17905159  2.29307897  0.11814451 -0.4371595 -0.52291238 -0.05198191  0.81765945  0.47613328",
  "6  -0.47624689  0.05974994  1.90236182  1.54758106  0.13403865  0.3311792  0.68374552 -1.75323736 -0.20979317  1.37857014",
  "7  -0.78860284 -0.70459646 -0.10097489 -0.13315096  0.22101947 -2.0142105 -0.06082195  0.09932759  0.37816777  0.45623640",
  "8  -0.59461727 -0.71721816 -1.35984070 -1.75652740  1.64084617  0.2119804  0.63296071 -0.57185006 -0.94540883 -1.13558847",
  "9   1.65090747  0.88465050 -0.66476944 -0.38877986 -0.21905038  1.2366750  1.33551762 -0.97400958  0.85692301 -0.43564547",
  "10 -0.05402813 -1.01559258  0.48545998  0.08920722  0.16806538  2.0375740  0.00729009 -0.17990623 -0.46103834  0.34610362",
  "11  0.11924524  1.95529397 -0.37560287  0.84501300  1.16838387  1.3011760  1.01755864  1.01494317  2.41677335 -0.64704563",
  "12  0.24368743 -0.09031959 -0.56187636  0.96252797  1.05418102  0.7567748 -1.18843404 -1.99274849 -1.65104890 -2.15764634",
  "13  1.23247588  0.21453883 -0.34391723  0.68430943  1.14526311 -1.7267304 -0.72160444 -0.42727929 -0.46398724  0.88425082",
  "14 -0.51606383 -0.73852770  0.09049665 -1.39527435 -0.57746800 -0.6015067  1.51921771  0.11663728  0.82537986 -0.82947761",
  "15 -0.99250715 -0.57438869  1.59850877  0.84964305  2.00248273 -0.3520465  0.37738797 -0.89320757  0.51013255 -0.57356027",
  "16  1.67569693 -1.31701613 -0.08856511 -0.44655722  0.06670087  0.7035239 -2.05222282  0.33390294 -0.58948104  1.50390061",
  "17 -0.44116322 -0.18292539  1.08079950  0.17480270  1.86685184 -0.1056713 -1.36403745  0.41142992 -0.99678074 -0.77414493",
  "18 -0.72306597  0.41898240  0.63075412  0.07455118 -1.35090269 -1.2586486 -0.20078102 -0.03303616  0.14447570  0.84573154",
  "19 -1.23627312  0.32430434 -0.11363990  0.42816676  0.02098359  1.6844357  0.86577940 -2.46589819 -0.01430741 -1.26068288",
  "-1.28471572 -0.78153649 -1.53290200  0.02467498  1.24991457  0.9113913 -0.10188326  2.57145815 -1.79028124 -0.35454240"
  )

[Destix]

You do not need an extra package. Use pairwise.t.test() from the base R stats package:

data <- df %>% pivot_longer(cols = starts_with("Col"))
data$name <- data$name %>% factor(levels = c("Col1", "Col2", "Col3", "Col4", "Col5", "Col6", "Col7", "Col8", "Col9", "Col10"))

pairwise.t.test(data$value, data$name,  p.adjust.method="BH")

      Col1 Col2 Col3 Col4 Col5 Col6 Col7 Col8 Col9
Col2  0.82 -    -    -    -    -    -    -    -   
Col3  0.97 0.87 -    -    -    -    -    -    -   
Col4  0.82 0.97 0.82 -    -    -    -    -    -   
Col5  0.82 0.82 0.82 0.74 -    -    -    -    -   
Col6  0.74 0.82 0.77 0.82 0.66 -    -    -    -   
Col7  0.82 0.97 0.88 0.95 0.82 0.82 -    -    -   
Col8  0.82 0.90 0.82 0.97 0.74 0.82 0.89 -    -   
Col9  0.97 0.82 0.97 0.82 0.82 0.74 0.82 0.82 -   
Col10 0.82 0.95 0.82 0.97 0.74 0.82 0.95 0.97 0.82

P value adjustment method: BH 

[nviet]

I appreciate this response. But I don't think this provides the correction in the context of a partial correlation analysis. I may have worded this poorly with the list(xxxx) code. Let me clarify a bit. For example, if I wanted to run:

pcor.test(df$Col1, df$Col2, list(df$Col3, df$Col4), method = "pearson")
pcor.test(df$Col1, df$Col6, list(df$Col3, df$Col4), method = "pearson")
pcor.test(df$Col1, df$Col7, list(df$Col3, df$Col4), method = "pearson")

Then had an output below but three of them:

estimate p.value statistic n gp Method
1 -0.397631 0.02955589 -2.293146 32 2 pearson

Then taking those p_values and throwing them into an object and running lets say a p.adjust(object, method = "BH").

That would be independent of the analysis. Unlike running a:

corr.test(df, y = NULL, use = "pairwise",method="pearson",adjust="BH",
alpha=.05,ci=TRUE,minlength=5,normal=TRUE)

Which makes the correction relative to the context of the analysis.

[Destix]

I'm sorry but I'm still struggling to understand where you're going with this.

A pairwise Bonferroni-Holm test for correlation coefficients is what you get from
corr.test(df, y = NULL, use = "pairwise", method="pearson", adjust="BH")
which tests the correlation between x and y.

A partial correlation test like

pcor.test(x, y, z)

tests the correlation between x and y under the condition z.

Those are two completely different statistical models.

[nviet]

Thank you. I understand that. I want to use the (adjust="BH") like in corr.test but in pcor.test. Is there a package that would allow that?

The problem is:

Adjusted p-values depend on the context of the analysis, the data, and the statistical method being used. If you're performing correlation tests and using the corr.test() function, the adjusted p-values from that analysis:

corr.test(df, y = NULL, use = "pairwise", method="pearson", adjust="BH")

will differ than when you do the same without

corr.test(df, y = NULL, use = "pairwise", method="pearson")

manually taking those p-values, assigning them into a object like below

pvalue2 <- c(0.20, 0.50, 0.17, 0.32, 0.42)

p.adjust(pvalue2, method = "BH")

So with this problem in mind, and put into the context of pcor.test, which does not allow adjust="BH" within this function, the p-values generated by using:

pvaluepcor <- c(0.20, 0.50, 0.17, 0.32, 0.42)

p.adjust(pvaluepcor, method = "BH")

wont reflex the most correct correction because it is not account for the entire context of the partial correlation analysis like corr.test does.

[Destix]

I'm afraid you have some misunderstanding about the theory of p-value adjustment methods and when to use them. In short, you would only use them in a hypothesis test where your test statistic is not constant.

Of course, this is not the case for a partial correlation test. You would never do a p-value adjustment here. Hence, the pcor.test() function does not have the option to do so, and rightfully so.

In multiple comparison correlation tests, your test statistic is a function of the number of test pairs. In this case, the p-value has to be adjusted.