我知道有很多关于如何避免 R 循环的主题,但我无法理解如何矢量化我的迭代。
我有一个数据集,在这里我用 m 表示。我想用这个函数生成一个新矩阵,它将由数据 (m) 每一列的相关系数的 p.values 组成。
m<-matrix(rnorm(100),nrow=10,ncol=10)
sig.p<-function(x){
n= ncol(x)
p.values<-numeric(n)
p.values<-matrix(nrow=n,ncol=n)
for (i in 1:C){
for (t in 1:C){
p.values[t,i]<-cor.test(x[,i],x[,t])$p.value
}
}
p.values
}
sig.p(m)
我无法理解如何使用 mapply(如果是这样的话)。
任何人都可以提供有关如何矢量化这些迭代的建议(使用 mapply 或其他)
提前致谢!
塞萨尔
最佳答案
您可以使用 rcorr
来自 library(Hmisc)
library(Hmisc)
rcorr(m)$P
或使用
library(psych)
corr.test(as.data.frame(m))$p
或使用
outer
来自 base R
outer(1:ncol(m),1:ncol(m), FUN= Vectorize(function(x,y)
cor.test(m[,x], m[,y])$p.value))
基准
我尝试了一个较小的数据集(
100*100
)和一个稍大的数据集( 1e3*1e3
)。以下是功能: akrun <- function() {outer(1:ncol(m1),1:ncol(m1),
FUN= Vectorize(function(x,y) cor.test(m1[,x],
m1[,y])$p.value))}
akrun2 <- function(){rcorr(m1)$P}
agstudy <- function() {M <- expand.grid(seq_len(ncol(m1)),
seq_len(ncol(m1)))
mapply(function(x,y)cor.test(m1[,x], m1[,y])$p.value,M$Var1,M$Var2)}
vpipk <-function(){
n <- ncol(m1)
p.values<-matrix(nrow=n,ncol=n)
for (i in 1:(n-1)){
for (t in (i+1):n){
p.values[t,i]<-cor.test(m1[,i],m1[,t])$p.value
}
}
p.values
}
nrussell <- function(){
sapply(1:ncol(m1), function(z){
sapply(1:ncol(m1), function(x,Y=z){
cor.test(m1[,Y],m1[,x])$p.value
})
})
}
在
100*100
数据集 library(microbenchmark)
set.seed(25)
m1 <- matrix(rnorm(1e2*1e2),nrow=1e2,ncol=1e2)
microbenchmark(akrun(), akrun2(), agstudy(), vpipk(),
nrussell(), unit='relative', times=10L)
#Unit: relative
# expr min lq mean median uq max neval cld
# akrun() 257.2310 255.9766 252.2163 254.4946 248.9807 246.5429 10 c
# akrun2() 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 10 a
# agstudy() 255.5920 258.0813 253.5411 256.0581 250.4833 249.0503 10 c
# vpipk() 125.8218 126.3337 125.4592 126.8479 124.9835 124.1383 10 b
#nrussell() 257.9283 256.8480 252.5297 256.0160 250.8853 242.0896 10 c
如果我改变
1e2
至 1e3
(没有时间做 microbenchmark
,但这里是 system.time
system.time(akrun())
# user system elapsed
#403.563 0.751 404.198
system.time(akrun2())
# user system elapsed
# 3.110 0.008 3.117
system.time(agstudy())
# user system elapsed
#445.108 0.877 445.947
system.time(vpipk())
# user system elapsed
#155.597 0.224 155.760
system.time(nrussell())
# user system elapsed
#452.524 1.220 453.713
关于R:如何在此函数中避免 R 中的 2 个 'for' 循环,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/27709227/