cor_mat , cor.test , and adjusting p values

Question

I want to determine the correlation and p-value between one variable and all of the others in an R data frame, so I created a correlation matrix with cor_mat, and then cor_get_pval. The relevant p-values were all p<.0001. To test the accuracy, I used cor.test on my one variable and some of the others, individually. I got the same correlation coefficient, but very high p-values (p>.5).

I think the reason for this is that cor_mat doesn't adjust p-values for multiple hypothesis testing, and cor.test does, but then shouldn't the cor_mat p-values generally be higher than the cor.test values? In any event, to minimize the potential for Type I error (which is why p-values get adjusted), I increased the confidence level in cor_mat to 0.9999, and got the same very low p-values. Here is my code:

cor_crime2022_test <- cor_mat(Crime2022_Gini , method = "pearson" , 
                              conf.level = 0.9999)
pval <- cor_get_pval(cor_crime2022_test)

cor.test(Crime2022_Gini$TotIncidents , Crime2022_Gini$Gini2022)

Is my thinking about this correct? Should I rely on the p-values in cor_mat or cor.test? Should I be doing something else?

Answer 1

This is not correct. The two should give exactly similar results. Note that cor_mat calls cor_test which calls cor_test_xy that uses cor.test from stats package:

rstatix::cor_mat
function (data, ..., vars = NULL, method = "pearson", alternative = "two.sided", 
    conf.level = 0.95) 
{    ...
    cor_test(data, vars = vars, method = method, alternative = alternative, 
        conf.level = conf.level) %>% as_cor_mat()
}

---

rstatix::cor_test
function (data, ..., vars = NULL, vars2 = NULL, alternative = "two.sided", 
    method = "pearson", conf.level = 0.95, use = "pairwise.complete.obs") 
{
     expand.grid(y = vars2, x = vars, stringsAsFactors = FALSE) %>% 
        as.list() %>% purrr::pmap_dfr(cor_test_xy, data = data, 
        alternative = alternative, method = method, conf.level = conf.level, 
        use = use) %>% add_class("cor_test")
}

---

rstatix:::cor_test_xy
function (data, x, y, method = "pearson", use = "pairwise.complete.obs", 
    ...) 
{
    if (is_grouped_df(data)) {
        results % doo(cor_test_xy, x, y, method = method, 
            use = use, ...)
        return(results)
    }
    suppressWarnings(cor.test(data[[x]], data[[y]], method = method, 
        use = use, ...)) %>% as_tidy_cor() %>% add_column(var1 = x, 
        var2 = y, .before = "cor")
}

---

Here is an example to show that the two are equivalent using the trees dataset.

(a <- rstatix::cor_mat(trees, method = 'pearson'))
# A tibble: 3 × 4
  rowname Girth Height Volume
* <chr>   <dbl>  <dbl>  <dbl>
1 Girth    1      0.52   0.97
2 Height   0.52   1      0.6 
3 Volume   0.97   0.6    1  


attr(a, 'pvalue')
# A tibble: 3 × 4
  rowname    Girth   Height   Volume
  <chr>      <dbl>    <dbl>    <dbl>
1 Girth   0        0.00276  8.64e-19
2 Height  2.76e- 3 0        3.78e- 4
3 Volume  8.64e-19 0.000378 0   

Now using the cor.test:

cor.test(trees$Girth, trees$Height)

    Pearson's product-moment correlation

data:  trees$Girth and trees$Height
t = 3.2722, df = 29, p-value = 0.002758
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.2021327 0.7378538
sample estimates:
      cor 
0.5192801 

Take note that the correlation is equal and also the pvalue is the same as that given by cor_mat.

You can get the pvalue and estimate for all the combinations:

unname(trees)|>
   combn(2, \(x)list2DF(cor.test(x[[1]], x[[2]])[3:4]), simplify = FALSE)|>
   do.call(what = rbind)

       p.value  estimate
1 2.757815e-03 0.5192801
2 8.644334e-19 0.9671194
3 3.783823e-04 0.5982497