R中矩阵每列之间的KS检验

Question

我想知道是否有一个包或一种新的方法可以比使用2个循环更快地在矩阵的每列之间执行KS测试？谢谢

all <- c()
for (i in (1:(ncol(a)-1))){
for (j in ((i+1):ncol(a))){
res <- ks.test(i,j)
all <- rbind(all,res)
}
}

Answer 1

假设我们有一个5列矩阵，每个矩阵包含来自不同正态分布的10个平局：

set.seed(1)
a <- sapply(seq(1, 3, 0.5), function(x) rnorm(10, x))
a
#>            [,1]       [,2]      [,3]     [,4]     [,5]
#>  [1,] 0.3735462  3.0117812 2.9189774 3.858680 2.835476
#>  [2,] 1.1836433  1.8898432 2.7821363 2.397212 2.746638
#>  [3,] 0.1643714  0.8787594 2.0745650 2.887672 3.696963
#>  [4,] 2.5952808 -0.7146999 0.0106483 2.446195 3.556663
#>  [5,] 1.3295078  2.6249309 2.6198257 1.122940 2.311244
#>  [6,] 0.1795316  1.4550664 1.9438713 2.085005 2.292505
#>  [7,] 1.4874291  1.4838097 1.8442045 2.105710 3.364582
#>  [8,] 1.7383247  2.4438362 0.5292476 2.440687 3.768533
#>  [9,] 1.5757814  2.3212212 1.5218499 3.600025 2.887654
#> [10,] 0.6946116  2.0939013 2.4179416 3.263176 3.881108

我们可以使用combn获得2列的所有唯一组合，并将ks.test应用于这些列对，检索单个向量中的所有p值，如下所示：

p <- apply(combn(ncol(a), 2), 2, function(x) ks.test(a[,x[1]], a[, x[2]])$p.val)
p
#>  [1] 0.1678213427 0.0524475524 0.0020567668 0.0002165018 0.9944575548
#>  [6] 0.4175236528 0.0123406006 0.1678213427 0.0524475524 0.4175236528

为了更清楚地了解这些p值所指的比较，我们可以将结果存储在一个5 x 5矩阵中，该矩阵将与循环的结果相同。不过，请注意，我们只需要运行ks.test10次，而不是25次，因为对角线已知为p=1，并且矩阵是对称的：

m <- matrix(0, ncol(a), ncol(a))
m[lower.tri(m)] <- round(p, 3)
m <- t(m)
m[lower.tri(m)] <- round(p, 3)
diag(m) <- 1
m
#>       [,1]  [,2]  [,3]  [,4]  [,5]
#> [1,] 1.000 0.168 0.052 0.002 0.000
#> [2,] 0.168 1.000 0.994 0.418 0.012
#> [3,] 0.052 0.994 1.000 0.168 0.052
#> [4,] 0.002 0.418 0.168 1.000 0.418
#> [5,] 0.000 0.012 0.052 0.418 1.000

或者，您可以将比较列的名称提供给p：

names(p) <- apply(combn(ncol(a), 2), 2, paste, collapse = "-")
p
#>          1-2          1-3          1-4          1-5          2-3 
#> 0.1678213427 0.0524475524 0.0020567668 0.0002165018 0.9944575548 
#>          2-4          2-5          3-4          3-5          4-5 
#> 0.4175236528 0.0123406006 0.1678213427 0.0524475524 0.4175236528 

^{创建于2022-09-10，reprex v2.0.2}