4

任何人都可以帮助我加快一些代码：

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矩阵，例如：NM

    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）执行回归分析（OLS） mat，然后使用包中的adf.test函数tseries并存储p-value. 最终结果sym.pr是所有的对称矩阵p-values（但实际上它不是 100% 对称的，请参阅此处了解更多信息），但它就足够了。

使用上面的代码，在一个600x300矩阵（600 个观察和 300 个对象）上，大约需要 15 分钟。

我想可能只计算对称矩阵的上三角形，但不知道如何去做。

有任何想法吗？

谢谢。

4

1 回答 1

2

从一些虚拟数据开始

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
于 2014-02-04T09:42:40.820 回答