任何能够帮助我加速一些代码的人:
n = seq_len(ncol(mat)) # seq 1 to ncol(mat)
sym.pr<-outer(n,n,Vectorize(function(a,b) {
return(adf.test(LinReg(mat[,c(a,b)]),k=0,alternative="stationary")$p.value)
}))
哪里
mat
是 NxM
N
的矩阵观察和M
对象,例如: Obj1 Obj2 Obj3
1 . . .
2 . . .
3 . . .
LinReg
定义为:# Performs linear regression via OLS
LinReg=function(vals) {
# regression analysis
# force intercept c at y=0
regline<-lm(vals[,1]~as.matrix(vals[,2:ncol(vals)])+0)
# return spread (residuals)
return(as.matrix(regline$residuals))
}
基本上我正在对
Obj1, Obj2
中的每个对象组合(即 Obj2,Obj3
和 Obj1, Obj3
和 mat
)执行回归分析(OLS) ,然后使用 adf.test
来自 tseries
的函数包装和存放 p-value
.最终结果 sym.pr
是所有 p-values
的对称矩阵(但实际上它不是 100% 对称的,请参阅 here for more info ),不过它就足够了。使用上面的代码,在
600x300
上矩阵(600 个观测值和 300 个对象),大约需要 15 分钟..我想过可能只计算对称矩阵的上三角形,但不知道如何去做。
有任何想法吗?
谢谢。
最佳答案
从一些虚拟数据开始
mdf <- data.frame( x1 = rnorm(5), x2 = rnorm(5), x3 = rnorm(5) )
我会首先确定感兴趣的组合。因此,如果我理解正确,那么
mdf[c(i,j)]
的计算结果应该相同。和 mdf[c(j,i)]
.在这种情况下,您可以使用 combn
函数,以确定相关对。pairs <- as.data.frame( t( combn( colnames( mdf ),2 ) ) )
pairs
V1 V2
1 x1 x2
2 x1 x3
3 x2 x3
现在,您可以将函数逐行应用于这些对(为简单起见,此处使用 t.test):
pairs[["p.value"]] <- apply( pairs, 1, function( i ){
t.test( mdf[i] )[["p.value"]]
})
pairs
V1 V2 p.value
1 x1 x2 0.5943814
2 x1 x3 0.7833293
3 x2 x3 0.6760846
如果你仍然需要你的 p.values 回到(上三角)矩阵形式,你可以转换它们:
library(reshape2)
acast( pairs, V1 ~ V2 )
x2 x3
x1 0.5943814 0.7833293
x2 NA 0.6760846
关于R加速方阵的矢量化,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/21547421/