在使用 cor 函数测试 R 中的显着性后确定 P 值答案

【问题标题】：Determining P values after using cor function to test significance in R在使用 cor 函数测试 R 中的显着性后确定 P 值
【发布时间】：2019-03-26 18:58:31
【问题描述】：

我是 R 新手，正在使用 cor 函数（特别是 Spearman）来确定时间流逝（如第 1 列所示）与其他 100,001 列中变量频率的增加之间是否存在显着相关性在我的数据框中。换句话说，我正在测试第1列和第2-100,001列之间是否存在相关性。

因此，目前，我将 CSV 文件导入 R，并将其转换为数据集 (test_data_fix)。目前，我使用以下代码，它返回一个数据框，其中包含我所有的列标签以及所有相关的 Spearman 相关值：

x <- test_data_fix[1:100001] y <- test_data_fix[1] corrs_test <- round(cor(x, y, method = "spearman", use="complete.obs"), 3)

但是，我还试图找到 P 值并将它们作为我返回的数据框中的列之一。我知道，如果我使用 cor.test，我可以一次检查单个相关性的 p 值，或者使用 corr.test 返回每个可能相关性的值。但是，有没有一种方法可以只返回 p 值，以测试第 1 列和所有后续列之间的相关性。

【问题讨论】：

标签： r correlation

【解决方案1】：

您需要迭代。例如，此方法为您提供 p 值矩阵，类似于 cor 为您提供每个列-列组合的相关值。

myfunc <- function(i,j) mapply(function(a,b) cor.test(mtcars[[a]], mtcars[[b]])$p.value, i, j)
mt <- mtcars[1:5]
outer(seq_len(ncol(mt)), seq_len(ncol(mt)), myfunc)
#              [,1]         [,2]         [,3]         [,4]         [,5]
# [1,] 0.000000e+00 6.112687e-10 9.380327e-10 1.787835e-07 1.776240e-05
# [2,] 6.112687e-10 0.000000e+00 1.802838e-12 3.477861e-09 8.244636e-06
# [3,] 9.380327e-10 1.802838e-12 0.000000e+00 7.142679e-08 5.282022e-06
# [4,] 1.787835e-07 3.477861e-09 7.142679e-08 0.000000e+00 9.988772e-03
# [5,] 1.776240e-05 8.244636e-06 5.282022e-06 9.988772e-03 0.000000e+00

甚至更好，有名字（感谢@RyanD）：

outer(mt, mt, Vectorize(function(a, b) cor.test(a, b)$p.value)) 
#               mpg          cyl         disp           hp         drat
# mpg  0.000000e+00 6.112687e-10 9.380327e-10 1.787835e-07 1.776240e-05
# cyl  6.112687e-10 0.000000e+00 1.802838e-12 3.477861e-09 8.244636e-06
# disp 9.380327e-10 1.802838e-12 0.000000e+00 7.142679e-08 5.282022e-06
# hp   1.787835e-07 3.477861e-09 7.142679e-08 0.000000e+00 9.988772e-03
# drat 1.776240e-05 8.244636e-06 5.282022e-06 9.988772e-03 0.000000e+00

如果您只需要将一列与所有其他列进行比较，那么：

outer(1, seq_len(ncol(mt)), myfunc)
#      [,1]         [,2]         [,3]         [,4]        [,5]
# [1,]    0 6.112687e-10 9.380327e-10 1.787835e-07 1.77624e-05
outer(mt[1], mt, Vectorize(function(a, b) cor.test(a, b)$p.value)) 
#     mpg          cyl         disp           hp        drat
# mpg   0 6.112687e-10 9.380327e-10 1.787835e-07 1.77624e-05

...但是将其作为data.frame 中的一列应用是没有意义的：添加一列表明（例如）第一个返回的 p 值将与第一行中的其他值相关联，绝对不是这样的。

【讨论】：

谢谢你！我可以仔细检查此方法返回的 p 值与逐一计算的 p 值是否正确。例如，如果我执行“cor.test(1:32, mtcars$cyl)”（我被告知如何将第一列（此处为 mpg）与指定的“其他”列（在本例中为 cyl）进行比较）然后它返回一个 0.8753 的 p 值，这显然比上面返回的 6.112687e-10 高很多。还是我误解了什么？再次感谢您的宝贵时间！
cor.test(1:32, mtcars$cyl) 在 "cyl" 和 1 到 32 的文字向量之间运行测试，它如何表示 "mpg"？试试cor.test(mtcars$mpg, mtcars$cyl)。
不，个别测试的 p 值应该与这个矩阵匹配。

【解决方案2】：

如果没有数据，以下内容未经测试，但我相信它可以满足您的需求。

它使用sapply 将第 2 列到第 100001 列中的每一列作为 x 和第一列作为 y 运行测试。

cor_test_results <- sapply(test_data_fix[-1], function(x)
  cor.test(x, y = test_data_fix[[1]], method = "spearman")$p.value)

【讨论】：

【解决方案3】：

cor.test() 确实可以为您提供一个 p 值（尽管它可能与关系有问题）。

也就是说，在更大的层面上，考虑问问自己将 100,000 列分组到一个样本中是否有意义。完全有可能它们不仅代表不同的样本，而且还可能对不同的人群进行抽样（尽管很难说不知道数据）。

此外，如果您确实决定进行一对一比较，那么如果您不应用某种类型的多重测试调整（这会以您的检测能力为代价），那么您的结果将几乎无法解释真阳性）。

【讨论】：